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INTRODUCTION 

This  report  is  an  attempt  to  analyze  Nordsieck's  paper  [l]  with  the 
help  of  the  asymptotic  methods  developed  Toy  Henrici  in  [2], 

The  definitions  of  the  Adams -Bashforth  formulas  and  the  Adams -Moult on 
formulas  are  contained  in  Section  1. 

An  asymptotic  relation  for  the  errors  in  the  Adams  methods  is  derived 
in  Section  2.   Most  of  the  arguments  used  are  directly  suggested  by  [2]. 

Section  3  is  concerned  with  the  self-stabilizing  property  of  the 
Adams  methods . 

Section  h   shows  how  to  get  Nordsieck-like  formulas  from  multistep 
methods;  applying  this  to  the  Adams -Moult  on  method,,  one  gets  the  Nordsieck 
formulas;  a  modified  Nordsieck  formula  is  suggested „ 

Predictor  formulas  for  the  Nordsieck  method  and  the  modified  Nordsieck 
formula  are  examined  in  Section  5  and  interpreted  in  the  language  of  multistep 
methods . 

The  errors  involved  in  the  process  of  changing  the  step  are  analyzed 
in  Section  6. 

Section  7  deals  with  starting  procedures. 

Section  8  gives  some  conclusions . 

A  table  of  results  is  presented  in  Section  9° 

Appendix  1  is  a  complement  of  Section  2, 

In  Appendix  2,  a  formal  generalization  of  the  Nordsieck  method  is 
discussed;  it  turns  out  that  these  formulas  are  "almost  equivalent"  to  the 
classical  multistep  methods.   As  a  particular  case,  to  every  multistep  method 
corresponds  an  equivalent  method  using  less  "past  values . " 

No  serious  practical  computations  have  been  done  for  checking  the 
theoretical  results  of  this  report.   However,  some  verifications  of  the 
Nordsieck  method  and  of  the  modified  Nordsieck  method,  k  =  2,  with  the  dif- 
ferential equation  y1  =  y,  have  been  satisfactory. 


Note:   In  this  report,  all  variables  are  real. 
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SECTION  1.   THE  ADAMS  METHODS 
Let  y(x)  be  a  function  and  y'(x)  be  its  derivative. 

Adams -Bashforth  Formulas  (Predictor  Type) 

Let  P(x)  be  a  polynomial  of  degree  k  +  1  satisfying  the  relations 

P(xQ)  =  y(xQ) 

P'(x0  -  ih)  =  y'(x0  -  in),    i  =  0,1,  ...,k; 

th 
P(x  +  h)  is  the  result  of  the  k   Adams -Bashforth  formula.   It  is  easy  to  show 

that  there  exist  the  constants  (independent  of  h,  y,  x)  (3  ,  (3  ,  „..,  ft  such 

that 

P(xQ  +  h)  =  y(xQ)  +  h(p0y'(x0)  +  P^'CXq  -  h)  +  ...  eky'(xQ  -  kh)); 


if  y(x)  possesses  k  +  2  continuous  derivatives,  the  difference 
L(y(x  ),  h)  -  y(x  +  h)  -  P(xn  +  h)  satisfies  the  relation 

L(y(xQ),  h)  =  Dkhk+2y(k+2)(|  ),      xQ  -  kh  <  |  <  xQ  +  h, 


where  D,  is  a  constant . 
k 


Adams -Moult on  Formula  (Corrector  Type) 

Let  Q(x)  be  the  polynomial  of  degree  k  +  2  satisfying  the  relations 

Q(xQ)  =  y(xQ) 

Q*(xQ  -  ih)  =  y'(x0  -  ih),     i  =  -1,0,  ...,k; 

Q(x  +  h)  is  the  result  of  the  k   Adams -Moulton  formula  and  satisfies  a  relation 
of  the  form 
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Q(xQ  +  h)  =  y(xQ)  +  hCp^y'CxQ  +  h)  +  PQy'  (xQ  -  L)  ...  +  ^y'  (xQ  -  kh)) 


where  |3  ,  (3_.  ...,    (3  are  constant.   If  y(x)  possesses  k  +  3  continuous 
derivatives,  the  difference  L(y(x  )_,  h)  =  y(x  +  h)  -  Q(x  +  h)  is  given  by 

L(y(xQ),  h)  =  Ckhk+3(0,      xQ  -  kh  <  |  xQ  +  h, 


where  C,  is  a  constant 
k 
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SECTION  2. 
AN  ASYMPTOTIC  RELATION  FOR  THE  ERROR  OF  INTEGRATION  WITH  THE  ADAMS  METHODS 

Notations 

Any  letter  with  an  arrow  is  a  vector  of  dimension  m;  the  components 

m 
of  any  vector  y  are  denoted  by  y  ,  y ,    ...,    y  .    |y|   =  Z  |y. |  is  the  norm 

-p  "*  i=1 

of  y. 

We  shall  consider  the  integration  of  a  system  of  m  differential 
equations  on  the  interval  b5xSc>  however,  in  order  to  allow  the  possibility 
of  initial  values  for  x  <  b,  we  shall  give  the  assumptions  on  the  differential 
equations  for  the  interval  a  <  x  <  c,  where  a  <  b. 

h  denotes  the  step  of  integration  and  we  write  x  =  b  +  nh,  where  n 

n         ' 

is  a  positive  or  negative  integer.   For  any  function  y(x)  defined  on  the  complete 

interval  [a,c]  or  only  at  the  discrete  points  x  ,  we  use  the  notation  y(x  )  =  y  . 
'  *  n'  n     n 

If  t  is  a  variable,  0(t)  represents  a  quantity  such  that 

|0(t)|  <  w|t|.  where  w  is  independent  of  t;  t  may  be  restricted  to  some 

interval.   0(t)  is  a  vector  such  that  | |o(t)|   -  0(t).   We  shall  use  in  the 

following  notations  of  this  type: 

y(x)  -  z(x)  +  0(hP); 

we  always  suppose  that  0(h)  is  uniform  in  x  for  a  <  x  <  c,  i.e.,  there  exists 
a  constant  w,  independent  of  x  and  h  such  that  |0(br  )|  <  w|h  |  for  a  <  x  <  c. 

For  any  function  y(x),  y   (x)  denotes  its  i   derivative. 

Theorem 

Let 


|j  7<x)  -  ?(x,y)  (1) 


be  a  system  of  m  differential  equations  and 
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g±i(x,y)  =  — ,      i,J  =  1,2,..., m 

be  the  elements  of  the  matrix  G(x,y). 

— »     — > 
We  assume  that  for  any  y  and  y*,  i,j  =  1,2,  .  .  „,m,  a  <  x  <  c: 

a)  g. . (x,y)  exists  and  is  continuous  in  x  and  y; 

b)  |gi.  (x,y)|  <  -  ,      .  n  constant  f2) 

c)  I g, .(*,?*)  -  g  .(x,y)|  <if\\y*   -  y||,     if   constant  (3) 

Let  z(x)  be  a  particular  solution  of  (l);  we  assume  that  for 
a  <  x  <  c. 

d)  z(x)  possesses  (p  +  2 )  continuous  derivatives  (see  (5)); 


Let 


e)  g..(x,z(x))  has  a  continuous  derivative  in  x,  i,j  =  l;...;m 


L(y(x),h)  =  y(x  +  h)  -  y(x)  -  h(j3_iy'(x  +  h)  +  ...  +  ^y' (x  -  kh))  (k) 


be  any  formula  described  in  Section  1  for  k  >  1;  if  p   =  0,  it  is  an  Adams - 
Bashforth  formula;  if  3   /  0,  it  is  an  Adams -Moulton  formula;  let  C  and  p 
(p  >  2)  be  such  that 

L(y(x),h)  =  ChP+1y(P+l)(e  ),      x  -  kh  <  |  <  x  +  h         (5) 


Finally  let  u  satisfy  the  difference  equation 

u  ,.  =  u  +  h(p  ,f  ,,  +  p_f  +  ...  a?  .  )  (6) 

n+l    n    v  -1  n+1    0  n        k  n-k  v 

->    ->.   -»•  . 
where  f  =  f (x  ,  u  )  for  n  >  q:  q  is  a  given  integer  such  that  x  ,  >  a;  the 
n     x  n'  n        —  U7  ^  °  q-k  — 

initial  values  u  ,  f   ,,...,  f  ,  satisfy  the  relations 

q'   q-17     '   q-k 

f)  u  =z  +  h  6  +  0(h    ),  5  constant,  s  integer  >  1  (7) 

f   .  =  z1  .  +  h  T).  +  0(h    ).  i  =  1,2,  ...,k.  r\ .    constant,  r  integer  >  1, 
q-iq-i      l         "      '  '    '    '      i         '  — 

(8) 
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Remark  that  u  , .  u  _,  ....  u  .  are  not  necessarily  defined.   (6)  is 
q-1/   q-2'     '   q-k  J  v  y 

an  implicit  equation  for  u   ;  in  order  to  provide  the  existence  of  a  unique 

solution,  we  assume  (see  Appendix  l)  that  h  <  h,  where 

g)  s^ssnr^i ;  rsJi  ♦  ...  iekn  (9) 

Then,  under  hypotheses  a  to  g,  the  error  d(x)  =  u(x)  -  z(x)  is  given  by 


d(x)  =  ^(e^x)  +  0(h))  +  hS(en(x)  +  0(h))  +  hr+1(em(x)  +  0(h))   (10) 


where  e  ,  e  ,  e    are  the  solutions  of  the  systems  of  differential  equations 


^  ez(x)  =   G(x,z(x))eJ(x)  -  Cz(p+l)(x),   e^b)  =  0;         (ll) 


^   e n(x)  =  G(x,z(x))ei;[(x),  e  (b)  -  8;        (12) 


~  eXII(x)  =  G(x,z(x))eIlx(x),  ^III(b)  =  ^        ^13^ 

n  -  (^  +  p2  +  . . .  +  ek)\  +  Os  +  P3  +  . . .  +  PkH2  +  •  •  •  +  Pk\- 

Generalized  Theorem 

In  the  preceding  theorem,  we  have  assumed  that  z(x)  is  a  fixed  solution 
of  (l).   In  fact,  the  same  conclusions  hold  if  z(x)  depend  on  h,  provided  we 
add  the  two  following  assumptions: 

h)  there  exists  a  constant  U,  independent  of  x  and  h  such  that  for 
a  <  x  <  c,  0  <  h  <  h 


II?(P+2)(X)||  <U; 
i)  there  exists  a  constant  V  such  that 


5  (x,z(x))|  <  V,      i,j  =  1,  ...,m. 
-'-J 
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Remarks 

1)  The  assumptions  b  and  c  are  very  strong  since  we  ask  them  to  be 

->     — > 
valid  for  all  y  and  y*;  however,  it  is  possible  a  posteriori  to 

weaken  them,  since  in  the  proof  of  the  theorem  we  shall  use  only 

a  bounded  region  of  y. 

2 )  The  proof  of  the  generalized  theorem  is  almost  word  for  word  the 
same  as  the  proof  of  the  theorem  itself  and  will  consequently 
not  be  given  here. 

Proof  of  the  Theorem 

All  the  symbols  introduced  so  far  in  Section  2  keep  their  meaning 
in  the  following , 

Lemma  1 , 

-*     -> 
Under  assumptions  b  and  c,  for  a  <  x  <  c  and  any  y  and  y*,  we  have: 

a)  I |y*  +  y||  <  ||y*ll  +  I  l?l  I J 

b)  |  |f  (x,y*)  -  f(x,y>|  |  <  a\  |y*  -  y|  |; 

c)  f(x,y*)  -  f(x,y)  =  G(x,y*)(y*  -  y)  +  0(||y*  -  y||2) 

d)  I  |G(x,y)y*|  |  <  ft|  \j*\  \ 

Proof 

Lemma  la  is  a  direct  consequence  of  the  properties  of  the  absolute 
values .   For  each  component  f.(x,y)  of  f(x,y),  we  have  the  mean-value  theorem 

m 
f.(x,y*)  -  f.  (x,y)  =  Z  g,.(x,y**)(y*  -  y  )   for  some  y**;      (ik) 


n   m  n 

f.(x,y*)  -    f.(x,y)|  <-  Z  |y*  -  y.|  =  -  |  |y*  -  y||; 

j  —  -J- 


||f(x,y*)  -  f(x,y)||  <  a||y*  -  y||, 

which  proves  lemma  lb„ 
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In  order  to  prove  lemma  lc,  we  start  with  equations  (l^-);  by  (3): 

gi(j(x,y**)  =  gi;j(x,y)  +  o(||7**  -  711); 

:ince  y*  -  y.    =  0( | |y*  -   y| | )  and  | | y**  -y||<j|y*-y||,    it  follows  that 
J          J 


f   (x,7*)  -  f   (x,y)  =     E  g,.(x,y)(y*  -  y.)  +  0(||y*  -  711), 

which  is  the  desired  result  for  the  i   component. 

Lemma  Id  can  be  proved  in  the  same  way  as  lemma  lb. 

Lemma  2 

Let  y_k,  y_k+1.>  •••,  yQ,    y1}    J2,    •••  satisfy  the  inequalities 

y  .  <  y  +  h(b  ny  .  +  b_y  +  .  .  .  +  b.  y  ,  )  +  7   for   n  >  0:     (is) 
Jn+1  -  Jn    v  -l^n+l    0Jn         k^n-k'  -  *     v  ' ' 

0  <  y.t  <  Y;   i  =  1,2, ...,k;   yQ  >  0;   b.  >  0,   i  =  -1,2, ...,k;   7  >  0; 


let  B  =  b  ,  +  b^  +  ...  +  b.  and  suppose  that 
-10         k 


h  >  0,   B  >  0,   h<|g;  (16) 


then 


yn  <  (^(y0  +  hBY  +  7)  +  ^)e3nhB   for   n  >  0.  (17) 


Proof 


From  (l6),  it  follows  that 

1  +  hB  <  1.5; 
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Thx  '   1  +  r^hB  5  1  -  2hB  <   (1  +  hB)2;  (18) 

from   (l6)  and   (l8),    we  get 

o 

y     ,    <   (l  +  hB)    fy     +  h(b^y     +  b_,y     _    +    .  .  .    +  b.  y     ,  )   +  7} ; 
Jn+1  -   v  lJn  v    Crn  lrn-1  kJn-ky  J ' 


since  y    .    <  Y  for   i   =  1,2,  .  .  .  ,k,   we   can  write  for  n  =  0,1,..  . ,k-l: 


y      _    <    (1   +  hB)    {y     +    (b_y     +  bny  +    ...    +  b  y    )   +  hBY  +   7}? 

•'n+l  —  n        v    On  l^n-1  n  0  •" 


from  the  two  last   inequalities,    it  follows  by  induction  that 


y     <  v     for  n  >  0,   where  v.,   =  y~,    and 
Jn  -     n  -     '  0       J0' 


v     .    =    (1  +  hB)    O     +  h(b_v     +  b_v     _    +    ...   b  v.)   +  hBY  +  7}  for  0  <  n  <  k-1, 

n+1  n  0  n  1  n-1  n  0  —       - 


v     .    =    (l  +  hB)    fv     +  h(b_v     +  bnv     .    +    . . .   b  v     .  )   +  7}  for  n  >  k 

n+1  L    n  0  n  1  n-1  k  n-k  '  - 


it  follows  that  v     <  v     <  v      ...;    consequently 


v     ,    <  (l  +  hB)2f(l  +  hB)v     +  hBY  +  7}  for  0  <  n  <  k   -   1, 

n+1  —  n  —       —  ' 


<   (1  +  hB)2{(l  +  hB)v     +  7}  for  n  >  kj 


n+1  -  L  v  n 


2  ^ 

by  induction  again,    since    (l  +  hB)     <  3>    (l  +  hB)     <  ^-,   we  have 

v     <  w     for  n  >  0,   where  w.   =  y^  and 
n  -     n  -     y  0       J0 


w         =  4w     +  3(hBY  +  7)  for  0  <  n  <  k  -  1;  (19) 


w      ,    =   (1  +  hB)3w     +  37        for  n  >  k;  (20) 

n+1        v  /      n  -     ?  v 
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(19)  is  a  difference  equation  the  solution  of  which  is  given  (as  it  is  easy 
to  prove  by  induction)  by 


4n   1 
w  =  4  w„   +  3   i  "  n  (hBY  +7)      for  0  <  n  <  k: 
n      0      4-1  —   —  7 


in  particular,  for  n  =  k 


w  =  4kw  +  (4k  -  l)(hBY  +  7)  <  ^k(w  +  hBY  +  7);         (2l) 


in  the  same  way,  the  solution  of  (20 )  is 


,.    ^3(n-k)      (1  +  hB)3(n"k)-l 
=  (1  +  hB)  v     w  +  - ■ - 


w   -  (1  +  hBj^v    '  w,  +  ■* '- 37      for  n  >  k 

n 


(1  +  hB)3  ■  1 
clearly,  we  can  write  for  n  >  0: 


f 


/  f-,        v,^^3n      0-v  (1  +  hB)3n  -  1 
<  (1  +  hB)   w.  +37  - -z i 


n  -  v      '  k 


(1  +  hB)3  -  1 


from  the  inequalities 


(1  +  hB)3n  <   e3nhB,      n  >  0, 


<  1 


(1  +  hB)3  -  1  "  3hB 
and  by  (21 ),  we  get 


w  ■:  (4k(wn  +  hBY  +  7)  +  r^)e3nhB,       n  >  0 
n  —  v   v  0  hB      '        — 


which  is  the  desired  result 
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Lemma  3 

If  a  function  y(x)  possesses  a  second  continuous  derivative  on 


a  <  x  <  c,  then 


L(y(x),h)  =  0(h2) 


Proof 

First  consider  the  case  y(x)  =  x„   By  (5)  it  follows  that  L(x,h)  =  0, 
i.e., 

P_x  +  P0  +  ...  +  3k  =  i;  (22) 

now,  let  y(x)  be  any  function  satisfying  the  hypothesis;  since  y'(x)  is  absolutely 
continuous  on  a  <  x  <  c,  y'  (x  +  9  )  =  y'  (x )  +  0(h)  for  any  a  <  x  -  kh  <  0   <  x  +  h  <  c 

L(y(x),h)  =  y(x)  +  h(y'(x)  +  0(h))  -  y(x)  -  hy'txKp^  +  PQ  +  ...  +  \  +   0(L)); 

by  (22),  L(y(x),h)  =  0(h2),  q.e.d. 

Lemma  h 

If  e(x)  is  the  solution  of  the  systems  (ll)  or  (12),  under  assumptions 
d  and  e,  then  e(x)  possesses  a  continuous  second  derivative. 

Proof 

By  lemma  Id,  the  system  of  differential  equations  have  a  Lipschitz 
condition  which  insures,  under  assumptions  d  and  e  a  first  continuous  derivative 
for  e(x)  (see  for  example  [2]).   The  lemma  h   results  by  differentiating  the 
equations  (ll)  or  (12). 

Lemma  5 

Let  t(x,y)  be  a  vector  satisfying  the  Lipschitz  condition 

I  I  t(x,y*)  -  t(x,y)  I   <  Il|  |y*  -  y|  I,      a  <  x  <  c,  II  constant. 

Suppose  that  the  function  y(x)  satisfies  the  system 


■11- 


—  y(x)  -  t(x,y(x)),      a  <  x  <  c, 


and  that  v  satisfies  the  difference  equation 
n 

v  n  =  v  +  h(B  ,t  .  +  p.t  +  ...  +  at  .  )  +  7   .  (23) 

n+1  n    v  -1  n+1    0  n         k  n-k     n'  v  J  ' 


— >    — *.        — *   >  ,  . 

where  t  =  t(x  ,v  j  for  n  >  j;  j  is  a  given  fixed  integer  (positive  or  negative) 

such  that  x  ,  >  a:  let  w  =  v  -  y  and  L(y(x),h)  he  the  vector  of  components 
j-k  —  '             n    n    n      ww'  * 


L(y,(x),h).   If 


1 


1)  I |w  I  j  <  W,  n  =  j-1, j-2, . . .jj-kj  W  constant 

2)  I 7      I    <  T,  n>j.r  constant 
1  '    n       —  —  w* 

3)  I  |L(y(x),h)|  I  A,      a  <  x  <  c,   A  constant 


then 


||w(x)||  <  (+k(||w.  II  +  hBw  +  r  +a)  +  L_±A)e3(c-a)B  for  x.  <x  <  c 
-      j  nB  j  -   - 

where  B  =  Il(  1(3^1  +  |pQ|  +  ...  IpJ)  and  h  <  |^. 


Proof 

If  we  subtract  the  relation 


y     ,    =  y     +  h(p  ny'    ,    +  B^y'    +   ...   +  B,  y'   ,  )  +  L(y  ,h) 

Jn+1        Jn  v    -lJn+l  CTn  KkJn-ky  wn'    y 


from   (23),   we  get 

w     _    =  w     +  h(p  _  (t     ,    -  y'      )  +   . . .   +  6,  (t     .    -   y'    ,  )   +  7     -   L(y  ,h); 
n+1  n  v    -lv    n+1        Jn+1  kv    n-k       •'n-k  n  wn'      f 

Since        t      -   y'        <  n|    w      |,   we  have  by  lemma  la: 
1  '    n        Jn       -      '  '    n'  " 

I  |w    .  1 1  <  I  |w  1 1  +  h(n|B  J   I  |w    .  1 1  +  ...  +  n|p.  I   I  |w    v|  I )  +  r  +  A; 

'      n+1       —  n1 '  '    — 1 '     ' '    n+11 '  ■   k'    '      n-k' '  ' 
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by  lemma  2,    if  we   set  y^  =    |  \^n+.  |  |,    7   =  T  +  A,    Y  =  V,    B  =  Jld^J    +    ...     \\\), 
we  get  for  n  >  0 


I -*        ll^/),k/||-*|i        1T1T7       _        A  \        T  +  Av   3nhB 
w      .        <   (4   (     w.        +  hBW  +  r  +  A )   +  — -=— )e  : 

n+j ' '    -  ' '    j  hB  ' 


since  a<x.<x      .<c,    it  follows  that  nh  <  c   -   a,    q.e.d. 
J  -     n+j  -  - 

Proof  of  the  Theorem 

We  first  consider  a  particular  case.   Suppose  that  for  i  =  1,2,,,,  „.<,k 

.   -> 
1}  u   .  are  defined; 
q- 1 

2)  f   .  =  f  (x   .oU   .  ); 
q-i      q-i   q-i 


3)  u   .  =  z   .  +  0(h  ):    (see  7). 
'        q-i     q-i        "    K 


From  the  third  condition,  because  of  lemma  lb,  we  have 


f   .  =  z'  .  +  0(hb). 
q-i    q-i 

Since  by  lemma  lb,  f(x,y;  satisfies  a  Lipschitz  condition  with  constant  fl,  we 

— >    — >         — >.    .    — >.    .   — »     ->   — >     —> 
can  apply  lemma  5^  replacing  t  by  f ,  II  by  Q,}    y{x)   by  z(xj,  v  by  u  ,  w  by  d  , 

J  by  q,  r  by  0,  W  and  |  |w.  |  |  by  0(hS),  A  by  0(hP+  )  (see  (5));  we  get 

J 

||d(x)||  •:  (Uk(0(hs)  +  hB0(hs)  +  0(hp+1))  +  0(hP^l)e3(c-a)B  =  Q(hs)   +  o(hp)  = 

—  LID 

(2k) 
Now,  let  us  subtract  from  (6)  the  corresponding  equation  for  z(x): 

d  _  =  d  +  h(p  . (f  .  -  z1  . )  +  p.(f  -  z' )+...+  ft  (f  .  -  z'   ))  -  L(z  ,h); 
n+1    n      -lv  n+1    n+1     0     n    n  k  n-k    n-k      v  n'  " 

by  lemma  lc,  for  n  >  q  -  k: 

f   -  z1  =  f(x,u  )  -  f(x,z  )  =  G(x  ,z  )d  +  0(   d    ). 
n    n      '  n       '  n       n'nn    Ul  n" 
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We  introduce  the  notations 

G(x)  =  G(x,z(x));  Gn  =  G(xn). 

By  (2^)  we  can  write  for  n  >  q 

d  .  -  d  +  h(p  ,G  nd  _+p.G  d  +  ...  +B  G  .d  ,  )  -  f(z  ,h)  +  0(h2p+1)  +  0{h2s+1); 
n+1    n      -1  n+1  n+1  0  n  n        k  n-k  n-k     v  rr       v      '     v     '> 

on  another  side  by  assumption  d,  (5)  and  (22),  we  have 

L(z  ,h)  -  ChP+1(B  nz(P;l}  +  ...  +  p  z(p+l})  +  0(hP+2); 
x  n'  -1  n+1  n  n-k  " 


then 


d  ,  -  t    +  h(p  .(g  nd   -  ch^P;^  +  ...  a  (G  vd  v  -  chpz(p+l}) 

n+1  n  v    -lv    n+1  n+1  n+1  Kk^    n-k  n-k  n-k 


P+2\    ,    7t^2p+lx    ,    ?t/v,2s+l- 


+  O(h^)  +  0(h^x)  +  Oth"^). 


We  have  to  distinguish  three  cases : 

->     -s— *       ~~* 
First  case:   p  >  s.   Let  d*  =  h  d  :  d*  satisfies  the  difference  equation 
n       n7      n 

d*  .  =  d*  +  h(p  nG  nd*+  ...  +p,G  .d*  .  )  +  0(h2), 
n+1    n      -1  n+1  n         k  n-k  n-k        " 


with  the  initial  conditions 

d*  =  6  +  0(h),  d*_  .  =  0(1),      i  =  1,2,  .  .  ,,k; 

we  compare  d*(x)  with  the  solution  e(x)  of  the  system  of  differential  equations 


f^  e(x)  =  G(x)e(x),      e(b)  =  5; 


Let  w  =  d*  -  e  :  we  have 
n    n    n' 


-ik- 


w  =  0(h),  w .  =  0(1),      i  =  1,2,... ,k; 

by  lemmas  3  and  k,    L(e(x),h)  =  0(h  );  we  can  apply  lemma  5  with  t(x,y)  =  G(x)y, 
n  =  fi  (lemma  Id),  vn  =  d*,  yn  =  en,  W  =  0(l),  r  =  0(h2),  A  =  0(h2),  j  =  q, 
w    =  0(h):  one  gets 


|  |w(x)|  |  <  (^(0(h)  +  hBO(l)  +  0(h2)  +  0(h2)  +  0(h2)h;°(h2)-)e3(c-a)B  =  0(h), 

i.e.,  d*(x)  =  e(x)  +  0(h)  and  d(x)  =  hS(e(x)  +  0(h)),  x  <  x  <  c. 

->■  -p-»   -> 
Second  case:   p  =  s.   Let  d*  =  h  d  ,°  d*  satisfies  the  difference  equation 
n       n^   n 

d*   =  d*  +  h(p  n(G  >  -  Cz(p!l})  +  ...  +  MG  J*  .     -   Cz(p:l}))  +  0(h2); 
n+1    n    v  -lv  n+1  n     n+1  Kkv  n-k  n-k     n-k        \      '* 


using  the  same  arguments  as  in  the  first  case,  we  get 

d(x)  =  hP(e(x)  +  0(h)),     x  <  x  <  c; 


where  e(x)  satisfies  the  differential  equation 


^   e(x)  =  G(x)e(x)  -  Cz(p+l)(x);      e(b)  -  5. 


Third  case :   p  <  s.   This  case  is  identical  to  the  second  one,  except  that  we 
have  to  replace  5  "by  0.   Clearly,  the  following  relation  is  valid  for  the  three 
cases 

d(x)  =  hP(ex(x)  +  0(h))  +  hS(en(x)  +  0(h)), 

where  e  (x)  and  e   (x)  satisfy  respectively  (ll)  and  (12);  this  finishes  the 
proof  of  the  theorem  for  the  particular  case. 

Now,  it  is  easy  to  prove  the  theorem  for  the  general  case.   Subtracting 
from  (6)  the  corresponding  equation  for  z(x)  for  n  =  q  we  get 
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Vl   -   dq  +   h(p-l<fq+l    "    ZQ+l'    +   P0<fq    "    *q>    +    •••    +   Mfq-k    "    Zq-k»    "    L\>^> 

(25) 

by  lemma  lb,       |f  -    z!       | |    <  fi| |d        |  |  ;   by  lemma  la  and  equations    (5),    (l), 

(8),   we  have 

||d       1 1   <  o(hs)  +  h(|p  Jn||d     ,  ||   +  0(hS)  +  0(h2))  +  0(hP+1): 

by  (9),  1  -  h|p  .  |ft  >  .5;  consequently 

l|dq+1M  =  0(hS)  +  0(hr+1)  +  0(hP+1); 
using  this  relation  in  (25)  we  get  by  (7)  and  (8): 

Vi  =  h"*  +  *f*\  +  *4**%  +  •••  +  Vr+1\  +  3(hS+1)  +  3(h"2)  +  °(hP+1)- 

— * 

In  the  same  way,  we  can  obtain  for  d    : 

v =  v +  e2hr+1^i +  h^X +  ••• +  vr+1\-i +  3(hS+1)  +  °<hr+2> +  °(hP+1» 

=  hsf+hr+1(p1+e2)n1+...hr+1(pk_1+pk«k_1+hI'+1pkrlk+o(h3+1)+o(hr+2)+o(hp+1); 


repeating  the  same  process  for  d   }    ...,    d   ;  we  finally  get 


d       =  h  6  +  h     0   +  p   +  ...  +  p^  +  h     o   +  ...  +  pk)n2  +  ...  +  n     pki} 


7*/,s+lN  -^/.  r+2  v        -*/.p+ls 

+  0(h        )  +  0(h        )  +  0(h^      ) 

,  s^       tT+1-*  7?/i_s+l\        7*/,  r+l\        t?/vP+1n 

=  h  6  +  h       ^  +  0(h        )   +  0(h        )   +  0(h^      ), 
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dq+k_.  =  0(hS)  +  0(hr+1)  +  0(hP+1),      i  =  1,2,... ,k, 


We  are  now  in  position  to  apply  the  theorem  proved  above  for  the  particular 
case.   It  suffices  to  replace  q  by  q+k  and  to  remark  that  the  system  of 
differential  equations 


^eIX(x)  =  G(x)en(x) 


is  linear, 
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SECTION  3-   SELF- STABILIZING  PROPERTY  OF  THE  ADAMS  METHODS 

We  use  the  same  notations  and  conventions  as  in  Section  2.   Let 
y  =  f (x,y)  be  a  system  of  m  differential  equations  satisfying  the  Lipschitz 
condition 

| |f(x,y*)  -  f(x,y)| |  <  a\  |y*  -  y| |      f or  a  <  x  <  c.        (l) 

->.  .  -> 

Let  z(x)  be  a  particular  solution  of  the  system  and  u  the  approximation 

obtained  by  using  the  k   Adams-Moulton  formula  (see  Section  l): 

u  .  =  u  +  h(p  ,?  ,  +  . ..  +  a  f  ,  ).  (2) 

n+1    n      -1  n+1         k  n-k  '  v 


where  x  =  b  +  nh,  f  »  ?(x„,u^)  for  n  >  0,  0  <  h  <  nCj(  \a      \ — —    |Q  '  r\  1  H. 


:.u  ;iorn^u,  u^-n^  ._,  /  i  • 1 1  _  '  i  > 


We  suppose  that  the  initial  conditions  for  u  in  b  are  such  that 

u(x)  =  z(x)  +  0(h),   b  <  x  <  c;  (2a) 

let  v(x)  be  the  exact  solution  of  the  differential  equations  such  that 

v(c)  =  u(c);  v(x)  depends  on  h;  we  suppose  that  v(x)  possesses  k+3  continuous 

derivatives  and  that  there  exists  a  constant  W  such  that 

| |v^k+3^(x)| |  <  W      for  b  <  x  <  c,  h  <  h. 


Theorem 

For  any  fixed  integer  r  independent  of  h 

u(c  -  ih)  -  v(c  -  ih)  =  0(hk+3),      i  -  0,1,  .  ..,r        (3) 


Remark 

The  self-stabilizing  property  allows  us  to  apply  the  so-called 
"Milne's  method  for  evaluating  the  local  discretization  error"  even  if  the 
initial  conditions  are  "bad." 
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Let  x  =  c  (t  depends  on  L)  and  u*  be  the  value  of  the  predictor 

\th 
obtained  by  the  (k+l)   Adams -Bashforth  formula 

H  "  Vi  +  h(p5?t-i  +  pI?t-2  +  •  •  •  +  ^.m'*  (*) 

in  Section  1,  we  have  mentioned  that 

vt  =  vi +  «tfu ♦  -  ♦  pfca-i.M) ♦  wk+37(k+3)<*>     (5) 

Subtracting  (5)  from  (k),   we  get  by  (l)  and  (3): 

in  the  same  way,  comparing  (2)  with  the  analogous  equation  for  v(x): 

vvvrvrc/3'(kt3)(^5^'  <?> 

where  x .  ,  ,  <  £  ,  f]  <  x ,  . 

t-k-1  —  ~  '     —  t 

Subtracting  (7)  from  (6): 

u*  -  vt  -  h(k+3)(Ckv(k+3)(,)  -  Dk+1v(k+3)(5))  +  0(hk+U);      (8) 

if  we  suppose  that  v     (x)  is  continuous  in  x  and  bounded  on  b  <  x  <  c, 
0  <  h  <  h,  then 

ck7(k+3)(^)  -  \+/k+3)(i)  -   (ck  -  Dk+1)7(k+3)(,)  +  3{h); 

replacing  in  (8),  one  gets  an  approximative  value  for  v     (Tl)  which  provides 
an  asymptotical  evaluation  of  the  truncation  error  C  h   v     (n)  in  (7)« 
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Proof  of  the  Theorem 
Lemma 

Under  hypotheses  (l)  and  (2a), 

u(x)  -  v(x)  =  0(h),   b  <  x  <  c. 

Proof 

By  (2a),  it  suffices  to  show  that  z(x)  -  v(x)  =  0(h).   By  (l),  we  can 
write 

x 

z(x)  -  v(x)  =  z(c)  -  v(c)  +  /  (f(z(u),u)  -  f(v(u),u))du; 

Jc 

x 

I  |z(x)  -  v(x)|  I  <  ||z(c)  -  v(c)||  +  /  I  |f(z(u),u)  -  f(v(u),u)||du 

Jc 

c 

<  I  |z(c)  -  v(c) I  I  +  n /   |z(u)  -  v(u) I |du; 

x 

from  this  last  inequality,  it  is  known  that 

||z(x)  -  v(x)||  <  ||z(c)  -  v(c)||efi(c"x). 

Since  z(c)  -  v(c)  =  z(c)  -  u(c)  =  0(h),  the  conclusion  follows. 

Now  we  go  to  the  proof  of  the  theorem  itself.   Subtracting  the  two 
equations 

VVith(p.iV  •••  *  Vt-i-i''  (9) 

\  "  Vi  +  h^-A  +  •  •  •  +  \%,*i)  +  ckhk+37(k+3)(i ),  (10) 

where  x  =  c,  x      <  f  <  x  ,  we  get  by  the  lemma  and  (l) 
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0  =  ut_x  -  ?    +  0(h2)  +0(hk+3) 


i.e. 


*  i^    1   i   Af  V>     —' 


VX- Vi   =  °(h  >  +°(*   }"  (11) 

Changing  t  in  t-1  in  (9)  and  (10),  one  gets  because  of  (ll) 


ut_2  -  v  '   =  0(h2)  +  0(hk+3); 


"by  induction,  it  follows  that  u    -  v    =  0(h  )  +  0(h    )  for  i  =  0,1,  ...,s, 

Li  *~  X       0  —  X 

where  s  is  any  fixed  number  independent  of  h;  with  this  result  in  mind,  we 
subtract  again  (10)  from  (9)  and  get  this  time: 


3>  ,  T*fJs*3 


0  -  u    -  v    +  O(h^)  +  O(h^); 


by  induction  it  follows  that 

k+3 


u, 


.  -  v    =  0(hJ)  +  0(h  °)      for  i  -  0, 1,2, . . ., s-k-1 

X       Xi  ~*  X 


Again  by  induction  for  j  =  2}J,,k}  .  .  . 

=  0(hJ)  +  0(hk+3)      for  i  -  0,1,2,...,  s-('j-2)(k+l); 


u,  .  -  v,  . 
t-i    t-i 


for  J  =  k+3,  we  get  (3) 
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SECTION  k.      THE  NORDSIECK  METHOD 

Note 

In  order  to  simplify  the  notations  and  avoid  confusion  with  different 
kinds  of  vectors,  we  shall  consider  in  this  section  and  in  the  following  ones 
only  scalar  differential  equations;  generalizations  for  systems  are  always  easy 
and  in  most  of  the  cases  obvious. 

Let  y1  =  f(x,y)  be  a  differential  equation  to  be  integrated  for  x  >  b 
with  the  help  of  the  general  multistep  method 


W  ■  Vn  +  °lVl  +   •••   +  a/n-j   +  h(p-lfn+l  +  Vn  +   •■•   +  Pkfn-k>>   n  2  ° 

(1) 


where  0c.}    $.}    k,  j,    (j  >  0)  are  the  constants  which  characterize  the  formula; 

some  of  the  a.'s  or  B.'s  may  vanish;  furthermore 
1       i  ' 


x  =  b  +  nh;  f  =  f(x  ,y  )      for  n  >  1; 

n         '   n      n'  n  —  ' 


(2) 


Y0,    y±,    >-,    y_y    f0,    t.mlf    •  -.,    f_k  are  the   initial  values;    fQ,    f_1, 


may 


be  different  from  f(xQ,y0),  f(x_1,y_1) 

We  operate  some  formal  transformation  on  (l).   Let  u  be  the  column 

-*■ 
vector  of  elements  y,y  ,,  .,.,   y   .,  f.f  n.  ....  f  ,  and  p  the  constant 

n'  •'n-l        n-j'   n'   n-1     '   n-k 

vector  of  elements  a_.  ol  ,  ....  a.,  hB_,  hB, ,  ....  hB,  .   (l)  can  be  written  in 
m       0'   1       J    0    1        k 


r-*T 


the  form  (p  being  the  transposed  of  p): 


y  .  =  p  u  +  hB  _f  . 
Jn+1   ^  n     -1  n+1 


u  .  =  Au  +  f  , a 
n+1     n    n+1 


fn+l  =  f(xn+ryn+l) 


y   for  n  >  0 


(3) 


M 


(5) 


where 
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'%    ai    •'•    aj|    hpo    hpi    • 

•■   hV 

10               0         1 

0 

0               10 

0 

0 

1  1 

0                |                  * 

1             0 

1       o . 

rvi 


a  = 


0 
1 
0 


the  matrix  has  been  subdivided  into  four  submatrices;  the  upper-left  submatrix 
is  square  of  order  j+1;  the  lower-right  submatrix  is  square  of  order  k+l„ 

For  given  u  ,  (3)>  (^);  (5)  define  an  iterative  process;  for  n  -  0,, 
(3)  and  (5)  allow  to  compute  f  and  (k)  u   ,    etc.... 

Let  T  be  a  regular  matrix  of  order  k+j+2,  the  elements  of  which  depend 
only  on  h;  we  consider  the  linear  transformation 


u  =  Tv  : 
n     n' 


replacing  in  (3),  (M,  (5)j  we  get  for  v  the  iterative  process 


-»T- 


y  .  =  q  v  +  hp  _f  . 
n+1      n     -1  n+1 


v   ,  =  Bv  +  f   _,b      >   for  n  >  0 
n+1     n    n+1       ' 


fn+l  =  f(Vyn+l] 


(6) 
(7) 
(8) 


■*   _T-*       -1.   -+    -l-> 
where  q  =  T  p,  B  =  T  AT,  b  =  T  a. 
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The  process  (6),  (7),  (8)  is  equivalent  to  the  process  (3);  (M;  (5) 
— >      — »  — >  — > 

in  the  following  sense:   if  u  and  v  are  chosen  such  that  u  =  Tv  ,  then  both 

processes  provide  the  same  values  for  y  and  furthermore  u  =  Tv  for  all  n  >  0, 

r  n  n     n  - 

In  order  to  get  a  particular  matrix  T,  let  us  consider  the  polynomial 
P  (x)  of  degree  k+j+1  satisfying  the  relations 


n 


P  (x   .  )  -  y   .      for  i  =  0,1,  .  ..,j; 
n  n-i'    n-i  '  7 


P'(x   .  )  =  f   .      for  i  =  0,1,  .  .  .,k; 
n  n-i     n-i  '  '    7    7. 


(9) 


this  polynomial  exists  and  is  uniquely  defined;  in  order  to  prove  that,  it 
suffices  to  show  that  the  only  polynomial  Q  (x)  of  degree  k+j+1  satisfying  the 
conditions 


S/Vi^  =  °'  i  =  °*1'-"'J 


Q*(x   .  )  =  0,  i  =  0,1,  ...,k 
tl   n-i      '       7    '         ' 


is  the  null  polynomial;  by  Rolle ' s  theorem,  there  exist  |  ,  |  ,  ...,  |.  such 
that 


x  >  |_  >  x  .  >  6_  ...  >  6.  >  x  ., 

n    1    n-1    2        j    n-j' 

Q^(ft±)  =  0,     i  =  i,...,j; 

QlCx)  is  of  degree  k+j  and  possesses  at  least  k+j+1  zeros;  consequently 
Q^(x)  s  0  and,  since  Q^xJ  =  0,  Q^x)  =  0,  q.e.d. 

Following  Nordsieck's  notations,  let  us  write  P  (x )  in  the  form 

'  n 

p    (    *          U)        h(    (2)^X   "  ^           (3)/X   "  XnV  (k+j+2)/X     '  Xn^+ 

P    (x)  =  vx        +  h(vv       I   — +  v  — )      +    ...    +  v 


nn  Kn\hjn\hy  n  V       h      J 

(10) 

"•  (l)       fk+i+2)  ~» 

let  v  denote  the  vector  of  components  v   .....  v       and  recall  that  u 
n  n  '    '  n  n 

has  the  components  y,  ...,  y   . ,  f  ,  ....  f  ,;  then  the  system  (9)  can  be  written 

n'     '   n-j'   n;     ;   n-k' 

in  the  form 

-24- 


Tv  =  u  : 
n    rr 


since  the  expressions  x   .  -  x  depend  only  on  h,  so  does  T^  T  is  regular,  since 
the  system  Tv  =  0  has  only  the  trivial  solution. 

The  particular  choice  of  T  we  have  got  with  the  help  of  the  polynomial 
P  (x )  gives  to  the  iterative  process  (6),    {l)}    (8)  two  advantageous  features: 

l)  Predictor  formula.   Because  of  (9)>  P  (x )  can  be  considered  as 
a  polynomial  of  approximation  of  the  true  solution  in  the 
neighborhood  of  x  ;  consequently  P  (x    )  provides  an  approximate 
value  for  y   ;  i.e.,  a  predictor  value: 

xW     \  (l)   ut    (2)    (3)  (k+j+lK 

P(x  ,)  =  vv+h(vv/+vw+...+vv  °         ); 
n  n+1     n        n      n  n       7 


2)  Modification  of  the  length  of  the  step.   Suppose  we  use  the 

step  h  for  x  <  x  and  the  step  h*  for  x  >  x  ;  all  the  quantities 

related  to  h*  will  be  denoted  with  a  star:   x*  =  x  , 

n    n' 

— *     — > 

x*   =  x  +  h*,  . . .;  q,  B,  b  also  depend  on  h  and  will  be 

be  replaced  by  q*,  B*,  b*;  it  is  natural  to  choose  v*  such  that 

P*(x)  =  P  (x),  where 
nv      nv  " 

n    n        x  n    V   h*  /    n    V   h*  J  " 


by  comparing  with    (10),    one  gets   the   simple   rule: 


(1)*  (1) 

v  =  v 


(2)*  (2) 

v  =  v 


(3)*       h*     (3) 

n  h      n     ' 


(k+j+2)*  ..   /h*\k+j         (k+j+2) 

n  Ui  n  ' 
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by  (9)  it  is  easy  to  interpret  the  meaning  of  this  rule  in 
terms  of  the  original  method  (l);  denoting  by  y*,  y*  ,  . .., 

y*  .,  f*.  ....  f*   the  initial  values  necessary  to  the 

n-j   n'    '   n-k  J 

integration  with  step  h*  from  x  ,  we  get 


y*  .  -  P  (x  -  ih*)      i  -  0,1, ..,,,}, 
n-i    nn      '  '  '    '    ' 


f*  .  -  P'(x  -  ih*)      i  -  0,1,  ...,k; 

n-i    n  n  '  '        7    7 


one  can  observe  that  generally  f*  .  ^  f (x   -  ih*. y*  . ) 

n-i    -  n        n-i 


The  Nordsieck  Formulas 

If  (l)  is  the  Adams-Moulton  method  described  in  Section  1  to  which 
one  applies  the  transformation  described  above,  then  (6),  (7),  (8)  represent  the 
Nordsieck  formula.   The  formula  corresponding  to  the  k   Adams-Moulton  method 
will  be  called  "k   Nordsieck  formula;"  the  corresponding  matrix  B  in  (7)  is 
of  order  k+2  and  the  corresponding  polynomial  P  (x)  which  is  of  degree  k+1  is 
defined  by  the  relations 

P  (x  )  =  y 
n  n     n 


P'(x  .  )  m  f      .,    i  =  0,1,  .  ..,k 
n  n-1     n-i'      '  7        ' 


The  Modified  Nordsieck  Formulas 

We  can  write  the  k   Adams-Moulton  formula  in  the  following  way: 


y  .  =  y  +  hpnf  +  hp  . f  .  , . .  +  hft  f  ,  +  Of  .  .  +  p  . f  . : 
n+1   ''n     0  n     -1  n-1        k  n-k     n-k-1    -1  n+1' 


we  simply  add  a  zero  term;  we  can  apply  the  same  transformation  on  this  formula 
as  we  did  for  the  Nordsieck  formula;  the  result  will  be  called  the  "k 
modified  Nordsieck  formula;"  the  corresponding  matrix  B  in  (7)  is  of  order  k+3 
and  P  (x)  which  is  of  order  k+2  is  determined  by  the  relations 
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P  (x  )  =  y 

nv  n'   Jn 


P»(x   .  )  =  f   .,  i  =  0,1, ...,k+l, 

n  n-i     n-i 


Remark 

The  k   Nordsieck  method  and  the  k   modified  Nordsieck  method  are 

both  equivalent  to  the  k   Adams-Moulton  method  in  the  following  sense:   denoting 

by  P  (x)  the  polynomial  corresponding  to  the  Nordsieck  method  and  P*(x)  the 

polynomial  corresponding  to  the  modified  Nordsieck  method,  then  all  three 

methods  will  give  the  same  values  for  y  ,  n  >  0  if 

n'   — 


P^x.i)  -  Pi*(*_±)  -  t_±, 


i  =  0,1,  .  ,.,k; 


the  value  of  PA*(x  ,  , )  has  no  influence  on  y  ,  n  >  0, 
0  x  -k-1'  Jn'        — 


Example 


For  k  =  1,  the  Adams-Moulton  method  is  given  by 


8h      _h        5h 
An+1   yn  +  12  n   12   n-1  +  12  n+lJ 


in  order  to  get  the  Nordsieck  formula,  we  have  the  sequence  of  intermediate 
stages : 


12'  ~  12      -1  "  12' 


-T 


u  =  (y  ,  f  ,  f  .); 
n    KJn'      n'      n-ly' 


A  - 


A 

8h 
12 

h 
12 

0 

0 

0 

\o 

1 
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12  ' 


a  = 


t>  f    \  d)    J  (2)/   "  n\  (3)/ 

P  (x)  =  v    +  h  iv     — +  v.  ' 


X  -  x  \     /^/x  -  x\2 
n  \     ( 

\   — Z —  +  v 
n       n         n   \   h    )  n 


-T    (    (1)    (2)    (3). 

n     n  '   n  '   n   7 

the  relations  P(x)=y.P'(x)=f.P'(x  .)  =  f  .  lead  to  the  relations 
n  n     n7   n  n     n7   n  n-1     n-1 


where 


T  = 


un  =  Tvn  (11) 


0 

°\ 

/l 

0 

0 

0 

1 

ol 

T-1- 

(° 

1 

0 

p 

1 

J 

Vo 

1 

2 

1 
"  2 

-T   -T_    ,.   7h  hs 
q   =  p  T  =  (1,  — ,  p 


Th  Ir 

12  Z 

'XAT  =  (   0     0  0   )  ,  b  =  T_1a 

-I  ° 


(6)  and  (7)  become  explicitly 


(l)  7h  (2)  h  (3)  5h  - 
V  =  v  +  - —  V  +  T  V  +  - —  f 
^n+l    n     12  n     6  n     12  n+17 


p(l)  =  y(l)  ^  7h  „(2)  ^  h  „(3)  ^  5h 

n+1    n 


12  n     6  n     12  n+1' 


n+1    n+17 


r(3)..iT(a)+i,      . 

n+1     2   n     2   n+1' 
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(1) 


.(2) 


(3) 


as  it  appears  also  from  (ll),  v^    -  y  and  v^  '   =   f  ;  replacing  uy* J   by  a  , 
one  gets  for  the  first  Nordsieck  formula 


7h  _    h      5h  „ 
yn+l  =  yn  +  12  fn  +  6  &n  +  12  fn+l' 


1  f    1  f 

an+l     2   n   2  n+1" 


(12) 


J 


In  the  same  way,  we  derive  the  modified  Nordsieck  formula  for  k  =  1; 
by  writing  the  Adams -Moult on  formula  in  the  form 


8h  _     h  _      __      5h  „ 
yn+l  =  yn  +  12  fn  -  12  fn-l  +  °fn-2  +  12  W 


we  get  successively 


-*T   f.      8h      h  nv       _  5h 


-KP 


=  (y  ,  f  ,  f   ,  f   ), 

n     n'   n   n-1'   n-2  ' 


A  = 


8h 
12 

h 
"12 

0 

0 

1 

0 

0 

1 

a  = 


the  polynomial  P  (x)  =  v'1'  +  h(v'2M — - — -J  +  v'3M — =- — -)     +   v 

n  n  vnVh/  n       \        h        /  nVh 


(U^L^V) 


satisfies   the   conditions   P   (x    )  =  y  .    P'(x      .)   =  f      ..    i  =   0.1,2;    it   follows 

n     n  n'      n     n-i  n-i'  '    '    y 

that 


u     =  Tv 
n  n 


(13) 


where 
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T  = 


1 

0 

0 

°\ 

0 

1 

0 

\ 

0 

0 

1 

-2 

3  / 

0 

1 

-k 

*/ 

T 


0 

0 

0 

1 

0 

0 

3 

1 

k 

-1 

h 

1 

1 

l 

Z 

"  3 

6 

-»T       -*Tm        /.,      7h       h         h- 


B  =  T     AT= 


q 

-  p  i 

-    K±j 

12' 

IS'    " 

"   k' 

i 

7h 

h 

£\ 

M 

12 

6 

+    \ 

\ 

/l2\ 

0 

0 

0 

1 

b  = 

--  t" 

l-» 

a  = 

" 

3 

l 

3 

1       }      J 

h 

2 

5 

/ 

\      *     / 

1 

1 

i  / 

r 

\    i/ 

6 

"  3 

2/ 

W 

(6)   and    (7)  become   explicitly 


(l)       7h      (2  h     (3) 

Jn+1  n  12     n  6     n 


5  Vn       +  12     n+1' 


(1)  (1)       7h     (2)       h     (3) 

n+1  n  12      n  on 


5  Vn        +  12     n+1' 


v^   =  f        , 

n+1  n+l; 


r(3)_       3v(2)       1     (3)  +  3v0O  +  3f 

T  _       —      -     T-     V  -      —     V  +T-V  +T-I  _,. 

n+1  4     n  2     n  4     n  k-     n+1' 


(M            1  (2)        1  (3)  1(h)        1    . 

n+1            6  n            3  n  2      n            6     n+1' 

setting  v^        =   a     and  v^  =  b    .    we  get  for  the   first  modified  Nordsieck  method: 
n              n              n              n' 


7h  „         h  h  ,  5h   _ 

y     .    =  y     +  —  f     +7^a     -rt     +f--f     , . 
Jn+1        Jn        12     n       6     n       +     n       12     n+1' 
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3  13  3 

a         =   -  f  f     -  ±  a     +rb     +  r  f     u 
n+1  4     n       2      n       4     n        4     n+1' 


b     ,    =  -  ?  f     -  -r  a     +  77  b     +zf        . 

n+1  on       3     n       2     n       6     n+1 


Another  Way  for  Deriving  the  Nordsieck  Formulas 

A  comparison  of  (12)  with  Reference  1  shows  that  Nordsieck  gets  its 
formula  in  a  slightly  different  form.   In  fact,  by  working  directly  on  the 
polynomial  P  (x),  one  obtains  the  "true"  Nordsieck  formula  much  faster. 

a)  Nordsieck  formula 

Let  p  <  q  be  integers  and  S(x:p,q)  be  the  Lagrangian 
polynomial  of  degree  (q  -  p): 


(x   -  x  )(x  -  x   ,  )  .  .  .  (x  -  x   ,  ) 

p/v P+l _ 3"1 

q 

th 


s(x:p'q)  =  (x   -  x')(x   -  x''")  •••  5    -  *  75  ' 

q    P   q    p+l        q    q-1 


the  polynomial  P  (x)  of  degree  k+1  for  the  k   Nordsieck  formula 
is  defined  by  the  relations 


p  (x  )  =  y  ,  (13) 


P'(x  )  =  f   .,    i  =  0,1,  ...,k;  (Ik) 

nv  n     n-i'        '  '    '  7 


P  ,(x)  will  be  determined  by 


P  (x  ,  )  =  y  , . 

nv  n+ly   Jn+1' 


P'   (x  .  .  )  =  f  _  .,    i  =  0,1,  ...,k;  (15) 

n+1  n+l-i     n+l-i'        ;  '    '  7  ' 


P'(x)  and  P'   (x),  polynomial  of  degree  k,  are  entirely 
determined  by  (ik)   and  (15);  one  checks  easily  that 
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P'  _(x)  =  P'(x)  +  (f  _  -  P'(x  .  ))s(x:n  +  1  -  k.  n  +  l):      (l6) 

n+1       n        n+1    n  n+1'   v  '  J'  v  ' 


let  fP  =  P*(x   . );  by  (10) 
n  n+1  ' 


fP  .  y(2)  +  2v(3)  +   (k)  (k  +  l)v(^). 

n       n       n  v       n    ' 


equalling  terms  of  some  degree  in  (l6)  allows  to  compute 

(2)        (k+2)    .  .        ..   ..     _   (2)        (k+2) 

v  ,,  ....  v'     as  linear  combination  of  v   .....  v    ' 

n+1'    '   n+1  n  '  '     n 

and  (f  _  -  fP). 
K  n+1 

It  remains  to  compute  v^  _.  =  P  .  (x   .  )  =  y  , ;  by 

n+1    n+lv  n+l'   "'n+l'  J 

Section  1,  y   .  =  Q(x    ),  where  Q(x)  is  the  polynomial  of 

degree  k+2  satisfying  the  relations 


Q(x  )  =  y  . 

n     n 


Q'(x  .)  =  f  .,      i  =  -1,0,  .  ..,k; 

n-i     n-i  '    '        '    * 


one  checks  easily  that 


Q'(x)  =  P'(x)  +  (f  .  -  P'(x   .  ))s(x:n  -  k,  n  +  l), 


Q(x)  =  Pn(x)  +  (fn+1  -  P;(xn+1))J  S(y:n  -  k,  n  +  l)dy; 

x 

n 


consequently 


-  „(1)  _ 


x  , 
n+1 


W  ■  v£i       Pn<xn+1>  +  (fn+l  "  fP)j  S(y:n  "  k'  n  +  l)dy      (17) 


X 

n 
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b)  Modified  Nordsieck  formula 

The  polynomials  P  (x )  and  P  , (x )  of  degree  k+2  for  the 
*     J  nv        n+lv        ° 

"th 
k   modified  method  are  determined  by  the  relations 


P  (x  )  =  y  ,           P'(x   .)  =  f   .,  i  =  0,1,  . ..,k+l; 

nx    n     n             n  n-i     n-i'  '  '    '    ' 

P  .  (x  ,)  =  y  ,,  P'   (x  _  .)  =  f  .  .,  i  =  0,1,  ...,k+l: 

n+1  n+1    Jn+1'  n+lv  n+l-i'    n+l-i'  '    '        '        ' 


then 


P;+1(x)  =  P;(x)  +  (fn+]_  -  fP)s(x:n  -  k,  n  +  l),  (l8) 


with  fP  =  P'(x  .)  =  v(2)  +  2v(3)  +  3^k)   +  ...  +  (k  +  2)v(k+3)i 
n  n+1     n       n       n  n    * 

equalling  the  coefficients  of  some  degree  in  (l8)  gives 

(2)        (k+3)    . .        ..   ..     _  (2)        (k+3) 

v   ',...,  v  -,    as  linear  combination  of  v   ,  ....  v      and 

n+17    '      n+1  n  '    '   n 

(f  .  -  fP). 
v  n+1 

By  Section  1,  y    =  Q(x    ),  where  Q(x)  is  the  polynomial 

of  degree  k+2  which  satisfies  the  relations: 


Q(x  )  =  y  ,      Q'(x   .)  =  f   .,      i  =  -1,0,.  ..,k; 
n    Jn  n-i     n-i'  '  '    '    ' 


then 


Q'(x)  =  P'(x)  +  (f  .  -  P'(x   .  ))s(x:n  -  k,  n  +  l); 


a  comparison  with  (l8)  shows  that  Q(x )  =  P  -.(x),  since 

Q(x  ,  )  =  y  ,  =  P  ,  (x  .  ) : 
v  n+17   Jn+1    n+lv  n+1 


P_-,(x)  ■■■■   Q(x)  -  Pjx)  +  (f_n   P:(x_J)J  S(x:n  -  k,  n  +  l)dx   (19) 

x 


Xn+1 
yn+l  =  FJ\^   +  <f„+l  -  fP)/s(x:n  -  k,  n  +  l)dx.  (20) 

X 

n 
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Remark 

A  conparison  of  (l6)  with  (l8)  and  (17 )  with  (19)  shows  that  the  k 

modified  Nordsieck  formula  is  identical  to  the  (k+l)   Nordsieck  formula,  except 

for  the  coefficient  of  (f   ,  -  f  )  in  the  expression  of  y  ..  which  is  the  same 
,  n+1  'n+l 

in  both  k   formulas . 


Example .  We  consider  again  the  case  k  =  1. 
Nordsieck  formula 


p  (x)  .  r*1'  ♦  h^'hr5  ♦  ^'{  __£   ); 

n       n       n   \   h   /    n   V   A   / 


(1)   .,  (2)/X  "  Xn\    (3) 

=  v    +  h(vv    : +  vA 

n        n   V   h  J  n 

by  (13)  and  (14),  v^  =  y  ,  v'2''  =  f  :  we  set  a  =  v^3';  by  (16) 


n'   n      n         n    n 


x-x  x-x  x-x 

P'   (x)  =  f  .  +  2a  .  r2i±  =  f  +  2a  -^—2  +  (f    -  fp)  — — » 

n+1       n+1     n+1     h       n     n   h        n+1     y   h 

it  follows  that  a  .  =  a  +  -(f  _  -  fP):  by  (17): 
n+1    n   2V  n+1     " 

x 

p  (y  -  x    )(y  -  x  ) 

y  .  =  P  (x  _  )  +  (f  .  -  fP)  /     ^~ —  dy, 

Jn+1    n^  n+17    v  n+1      'J  n12         J ' 

2h 
x 
n 


f 


y   ,  =  y  +  h(f  +  a  )  +  |£  (f   ,  -  fP); 
Jn+1   Jn    v  n    n'   12  v  n+1     " 


we  get  for  the  first  Nordsieck  formula  the  expressions 


y  _  =  y  +  h(f  +  a  +  t|  (f   ,  -  fP)), 
Jn+1   Jn    v  n    n   12  v  n+1     '" 


fP  =  f  +  2a  , 
n     n 


a    =  a  +  \   (f  .  -  fP) 

n+1    n   2  v  n+1 


•3+- 


Modified  Nordsieck  formula 

As  before,  we  substitute:  v    =  y  »  v^  '  =  f  ,  v    =  a  . 

n      rr   n      rr   n      n 

v^  '  =  b  .   By  (19): 

n      n 


/x-x  \    /x-x  \2    /x-x  \3  r  (y-x  .,  )(y-x  ) 

■  ^n ff  %  -r    *n(-Ts)  ^(^1-*p) J   ""Jg  P  *■ 

^  '  x 

n 

by  taking  the  successive  derivative  at  x    for  the  two  members^ 
one  gets 


where 


y    =  y  +h(f  +a  +b  +  ^  (f    -  fP)), 
n+1    n      n    n    n   12  v  n+1 


=  2a^  +  3b  +  2   (f  _  -  fP), 


n+1     n     n   4   n+1 


n+1    n   0   n+1      " 


f P  =  f  +  2a  +  3b  . 
n     n     n 
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SECTION  5.   PREDICTOR  FORMULAS 

a)  Nordsieck  formula 

For  the  k   Nordsieck  formula,  the  predictor  value  for  x  _,  is  given 

'  n+1    & 

by  P  (x    ),  where  P  (x)  is  the  polynomial  of  degree  k+1  determined  by  the 
relations 


P  (x  )  =  y  ,      P  (x   .)  =  f  ..      i  =  0,1,  ...,k; 
n  n     n'       n  n-i     n-i  *  '    ' 


a  conparison  with  Section  1  shows  that  P  (x    )  is  nothing  else  then  the  value 

th  n  k+2  (k+2) 

of  the  k   Adams -Bashforth  formula  with  a  truncation  error  D  h   y     (!)• 


b)  Modified  Nordsieck  formula 

th 
For  the  k   modified  Nordsieck  formula,  the  predictor  value  for 

x    is  given  by  P  (x   ),  where  P  (x)  is  the  polynomial  of  degree  k+2  determined 

by  the  relations 

P  (x  )  =  y  ,      P*(x   .)  =  f   .,      i  =  0,1,  .  ..,k  +  1; 

n  n     n'       n  n-i     n-i7  '  '    '■  7 

it  follows  that  P  (x   )  is  the  value  of  the  (k+l)   Adams -Bashforth  formula 

with  a  truncation  error  D  nh   y     (l)i 

k+l 


Remark 

The  modified  formula  has  the  advantage  that  the  predictor  and  the 
corrector  formulas  have  a  truncation  error  of  the  same  order. 
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SECTION  6.   MODIFICATION  OF  THE  LENGTH  OF  THE  STEP 

Suppose  we  integrate  the  differential  equation  y'  =  f(x,y)  from  a  to 
c  with  the  step  h  for  x  <  b  and  h*  for  x  >  b. 

Let  u(x)  be  the  discrete  approximative  solution  for  x  <  b  and  u*(x) 
be  the  approximative  solution  for  x  >  b.   Let  z (x )  be  the  exact  solution  of  the 
differential  equation  such  that  z(b)  =  u(b);  as  h  varies,  z(x)  varies  also. 

We  use  the  notations 


x  =  b  +  nh,       x*  =  b  +  nh*: 

n         '  n 


u  =  u(x  ),        u*  =  u*(x*) 
n      n  '  n     v  n  • 


f  =  f  (x  ,u  )      for  n  <  0: 
n      n'  n  — 


f*  =  f(x*  u*)      for  n  >  0: 

n    ^  n'  n  —  7 


z   =  z(x  ),        z*  -   z(x*) 
n      n  n      n 


We  suppose 

1)  p  =  — ■  is  a  constant  as  h  •+  0: 

h* 

2)  all  the  hypotheses  necessary  to  apply  the  self -stabilizing 
property  and  the  generalized  theorem  in  Section  2  are 
filled. 

Because  of  the  self-stabilizing  property,  it  is  suitable  to  compare 
u*(x)  with  z(x).   Let  d*(x)  =  u*(x)  -  z(x)  for  x  >  b.   Replacing  d(x)  by  d*(x), 
h  by  h*,  q  by  0,  the  generalized  theorem  gives  the  asymptotical  relation 


d*(x)  =  h^(ei(x)  +  0(h))  +  h*S(eIZ(x)   +  0(h))  +  hr+1(e1T1(x)   +  0(h)),   x  >  b; 


the  term  e  (x )  corresponds  to  the  accumulation  of  the  local  truncation  errors 
and  is  comsequently  independent  of  the  starting  procedure  for  the  new  h*;  the 
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term  e   (x )  corresponds  to  the  error  at  x  =  b;  by  hypothesis  d*(b)  =0,  i.e., 
s  =  _oo  and  the  second  term  of  the  asymptotical  relation  vanishes;  only  the 
term  e    (x )  depends  on  the  starting  procedure. 


Nordsieck  Method 

Consider  the  k   Nordsieck  formula;  according  to  Section  k,    P^(x)  is 
the  polynomial  of  degree  k  determined  by  the  relations 

P^(x_.)  =  t     ,  i  =  0,1,  ...,k; 


the  starting  values  f*.  for  the  new  h*  are  given  by 


f?i  =  Pc,(x?i)>      i  =  0,1,  ...,k; 


because  of  well-known  properties  of  the  interpolation  there  exists  the  constant 
coefficients  (independent  of  h*)  a. .  such  that 


f*.  =  E  a.  .f  ,,       i  =  0,1,. .,,k,  (l) 

-i   ,=0  iJ  -J 


z'*  =  Z  a. .z' .  +  R.,    i  =  0,1,  .  .  .,k,  (2) 

-1        ij  -J     x1 


where  the  interpolation  error  R.  is  given  by 

(*;,  -  *p )(*?!-  mj  •••  ^!j  -  *_k)    (k+2) 
\-- (k ;  i): z       bi.*'        i  =  o,i,...,k 

with  min(x*.,x   )  <  |   <  x  .   Since  h  =  ph*  and  z^k+2^(|  .  )  =  z^k+2'(b)  +  0(h), 
we  can  write 

R.  .  h*k+l  (-i)(p  -  i)(gp  ■  O  ...  (kp-  1)  z(k+2)(b)  +  0(h#k+2);    .  .  0jl;...jk. 

subtracting  (2)  from  (l),  one  gets 

-38- 


f*.  -  z1*  =  -R.  +  0(h*   ), 
-1-1     1    v      ' 


since  by  the  self-stabilizing  property  f  .  -  z' .  =  0(h*    ) 

J  J 


Using  the  notation 


.  (-i)(p  -  l)(2p  -  i)  ...  (kp  -  i) 
\  ~  (k  +  i): 


we  get  by  the  generalized  asymptotical  theorem 


d*(x)  =  h^k+2(ei(x)  +  0(h))  +  h*k+2(e;[II(x)  +  0(h)),   b  <  x  <  c     (3) 


where 


^  ex(x)  =  G(x,Z(x))eI(x)  -  Ckz(k+3)(x),   e^b)  =  0, 


dx  eHI('X^  =  G(x>z(x))eiii(x)> 


eIi;r(b)  =  -  {(px  +  P2  +  ...+  Pk)Ex  +  (02  +  P3  +...+  Pk)E2  +...+  PkEk)z(k+2)(b), 

th 
we  recall  that  |3  ,    f3  ,  P  ,  .  .  . ,  p  are  the  coefficients  of  the  k   Adams- 

~ _L     U     _l_  K. 

Moult on  method. 


Remarks 

1)  The  error  due  to  the  modification  of  the  step  is  non-negligible  in  comparison 
with  the  error  due  to  the  accumulation  of  the  truncation  errors. 

2)  In  (3),  G(x,  z(x))  is  a  scalar,"  for  systems  of  differential  equations,  it 
suffices  to  replace  d*(x),  e  (x),  e   (x),  z(x)  by  the  corresponding  vectors 
and  G(x,z(x))  by  a  matrix. 
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Modified  Nordsieck  Formula 

Consider  the  k   modified  formula;  according  to  Section  k,    PJL(x) 
is  a  polynomial  of  degree  k+1  determined  by  the  relations 


Pq(x_.)  =  f_±,  ±  =   0,1,.,,,  k+1; 


the  initial  values  f*.  for  the  new  step  h*  satisfy  the  relations 

f!i  =  F^x-^'  i  =   0,1,. .„k; 


there  exist  the  constant  coefficients  b . .  (independent  of  h*)  such  that 

k+1 


f*  =  Z  b,,f  ,,     i  =  o,i,,.,,k,  (k) 


j=o 


k+1 


z'*  =     Z  b.  .z'  .  +  r.,      i  =  0,1,  ...,k,  (5) 

-i     Q  ij  -j    i'  ' 


where  R.  is  the  interpolation  error 

„    (X-1  -  X0)(X-1  -  Xl}  •••  (X-1  -X-k-l}   (k+3)„  ,       ,   „  ,     v 

with  min(x*.,  x  .  , )  <  6 .  <  x_.   We  can  write 
-i'   -k-1  —  "  i  —  0 


R.  =  h*k+2H.z(k+3)(b)  +  0(h*k+3), 


II. 


-  (-i)(P  -  i)  ...  ((k  +  l)p  -  i) 


1  (k  +  2).' 

because  of  the  self-stabilizing  property,  if  we  subtract  (5)  from  (h)t    one  gets 

f*  -  z'*  =  -R.  +  0(h*h+3%- 
-l    -l     l    v      ' 
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by  the  generalized  asymptotical  theorem,  one  obtains  the  relation 


where 


d*(x)  =  h*k+2(ei(x)  +  0(h))  +  h*k+3(eni(x)  +  0(h))         (6) 


gj  ei(x)  =  G(x,z(x))  -  Ckz(k+3)(x),      ex(b)  =  0; 


dx  eHI^X^  =  G(x^z(x))eIII(x)^ 


eIT1(h)   =   -    {(P1  +  32  +  .  ..  +  Pk)H1  +  (P2  +  ...  +  Pk)H2  +  ...  +  p^Jz^^^b); 

th 
we  recall  that  8  _  ,  8_.  8.  .  ....  6n  are  the  coefficients  of  the  k   Adams-Moulton 

-1'  0'      1'       k 

formula . 


Remarks 

1)  Except  for  x  closed  to  b,  the  error  due  to  the  modification  of  the  step  is 
negligible  in  (6); 

2)  the  same  relation  (6)  holds  for  systems  of  differential  equations  if  one 
replaces  the  scalars  by  the  convenient  vectors  or  matrix. 
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SECTION  7-   STARTING  PROCEDURES 

Although  the  starting  procedure  is  not  the  most  characteristic  feature 
of  the  Nordsieck  method,  it  is  most  painful  to  analyze.   In  order  to  make  the 
job  easier,  we  shall  introduce  a  so-called  "simplified  starting  procedure." 

Given  a  differential  equation 

y'  =  f(x,y)  (1) 

and  initial  value 


the  problem  is  to  find  a  "good"  initial  v~  (or  equivalently  an  initial  poly- 
nomial P  (x) ) . 

We  introduce  the  notations : 


x  =  b  +  nh,    f  =  f(x  .y  ),    n  >  0; 

the  norm  | |x|   of  any  vector  x  is  the  sum  of  the  absolute  values  of  its  components 
We  suppose : 

1)  f(x,y)  satisfies  the  Lipschitz  condition  |f(x,y*)  -  f(x,y)|  <  ft|y*  -  yrj., 
b  <  x  <  c; 

2)  the  solution  z(x)  of  (l)  with  z(b)  =  y_,  possesses  (k+k)   continuous 
derivatives. 


Starting  Procedure  for  the  k   Nordsieck  Method 

Starting  with  any  polynomial  P  n(x)  such  that  P  n(xn)  =  yn  and 
P'  n(xn)  =  fn,  k  steps  of  integration  are  executed  forwards;  then  h  is  replaced 
by  -h  (by  the  standard  procedure  of  modification  of  the  step)  and  k  steps  are 
executed  backwards;  -h  is  replaced  by  h;  let  Q  n(x)  he  the  polynomial  obtained 
at  that  stage;  the  new  polynomial  P  ,(x)  is  determined  by  the  conditions 

\J  9    -L 
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po,i<*o>  =  V  pi,i(xo>  "  fo>  po,i(x)  ■  90,0<X'' 

the  same  procedure  is  repeated  with  P  -,(x);  one  gets  a  sequence  of  polynomials 

po,o(x)<  V(x)'  po,i(x)'  %,iM>  ■■'■■ 

th 
Simplified  Starting  Procedure  for  the  k   Nordsieck  Method 

Starting  with  any  polynomial  P  n(x)  such  that  Pn  n(xn)  =  yn  and 

P'  _(x_)  =  f_,  k  steps  of  integration  are  executed  forwards;  let  P.  ~(x)  "be 
0,0  0     Cr      ^         °  '      k,  0X 

the  polynomial  obtained  at  this  stage;  P  ,(x)  is  then  defined  by 

V<xo>  ■  V  V(x)  ■  pk,o(x)' 

the  same  procedure  is  repeated  for  the  new  initial  polynomial  Pn  , (x);  a  sequence 

po,o(x)'  Pk,o(x)>  po,i(x)>  Pk,i(x)>  •■•  is  defined' 

For  the  Nordsieck  formulas ,  we  have  the  following  theorem. 


Theorem  1 

For  h  sufficiently  small,  all  the  polynomials,  P_  . (x)  and  Q   .(x) 
for  the  starting  procedure,  P   . (x )  and  P   . (x )  for  the  simplified  starting 

U,  1  K,  1 

procedure,  coverge  as  i  goes  to  infinity  to  the  same  polynomial  U(x)  of  degree 
k+1  satisfying  the  relations 

U(xJ  =  yn,  (2) 


U'(x.)  =  f(x.,U(x.)),      i  =  0,1,.  ..,k.  (3) 


th 
Starting  Procedure  for  the  k   Modified  Nordsieck  Formula 

P_  . (x )  and  CL  . (x )  are  defined  in  the  same  way  as  for  the  Nordsieck 
0,i         0,i  J 

formula,  except  that  (k+l)  steps  of  integrations  are  executed  forwards  and 
backwards  at  each  iteration. 
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Simplified  Starting  Procedure  for  the  k   Modified  Nordsieck  Formula 

It  is  identical  to  the  simplified  starting  procedure  for  the  k 

Nordsieck  method  except  that  (k+l)  steps  of  integration  are  executed  at  each 

iteration  and  define  the  sequence  P.  ~(x),  P.  .  „(x),  P^  ,  (x),  Pn  .  n(x).  ... 

0,Cr  ''   k+l,Ov  "   0,lv  "      k+1,1 

For  the  modified  Nordsieck  formula,  we  have  the  following  theorem. 


Theorem  2 

For  h  small  enough,  all  the  polynomials,  P   . (x )  and  Q  . (x)  for  the 
starting  procedure,  P   . (x)  and  P    . (x)  for  the  simplified  starting  procedure, 
converge  as  i  goes  to  infinity  to  the  same  and  unique  polynomial  U(x)  of  degree 
k+2  satisfying  the  relations 

U(0  -  Yn> 


U'(x.)  =  f(x.,U(x.)),      i  =  0,l,...,k+l, 


Proof  of  the  Convergence  of  the  Starting  Procedures 

We  shall  consider  only  the  case  of  the  simplified  starting  procedure 
for  the  k   Nordsieck  method;  all  the  other  cases  can  be  treated  analogously. 

We  first  recall  that  v  is  the  vector,  the  components  of  which  are 
the  coefficients  of  the  polynomial 


p   ,    y  (1)        h,    (2){X    ~   Xn\  (k+2)^_l^n\k+\ 

P     (x)    =    Vv  +    h(V  : +     ...     +    V  : ) 

n  n  vn\h/  n  \h/ 


and  that  u  is  the  vector  of  the  components  y,f,f  ,.,..,  f  ,  .   Since 

n  *  Jn}      n}      n-1'     '   n-k 

P  (x  )  =  y  ,  P1 (x   .  )  =  f   . ,  i  =  0,1, . . . ,k,  there  exists  a  regular  matrix  T 
n  n    "'n'   nv  n-i     n-i'      '  ' 

such  that  u  =  Tv  .   Let  u   .  and  v  ,  be  the  components  of  P   . (x  ) 
n     n        n,  l      n,  k  n,  l 

(n  =  0, 1,  ...,k);  proving  the  convergence  of  P   .(x),  (or  P   .  (x))  for  i  -»  °° 

U .  1  K  «  X 

— >         ->  ' 

is  equivalent  to  the  proof  of  the  convergence  of  vn  .  or  un  . . 

—»  —*  .—*   . 

Let  for  any  u_,  k(u  ;  be  the  new  initial  vector  obtained  after  one 

iteration  of  the  simplified  starting  procedure.   From  Appendix  1,  it  follows 
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that  the  convergence  of  u_  .  towards  a  limit  vector  independent  of  the  initial 
vector  un  _  will  be  insured  if  we  can  show  that  there  exists  a  constant  K  <  1 
such  that  for  any  two  u  and  u*  (u!l     uq      Yq;  uq   =  uo      fo^  we  have 


k(u*)  -  k(uq)\  I  <  K||u*  -  uQ| 


We  apply  to  u  and  u*  the  simplified  starting  procedure;  the  quantities  con- 

— > 
cerning  u*  will  be  denoted  by  a  star.   Let 

y   ,  =  y  +  h(B  _f   ,  +  ...+  ft  f  .  ) 
Jn+1   Jn  v  -1  n+1         Kk  n-ky 


be  the  Adams -Moult on  formula;  for  y  and  y*,  we  get 

yi  =  yo  +  h^-ifi  +  pofo  +  "•  +  V-k*' 
yJ  =  vg  +  h(p_1f*+:P0fS+  ...  +Pkf?k); 

but    |f*  -   f    |    <  ft|y*  -   y    |;    if  1   -    |p      |hfl  >  0,    then 

lv*  -  y  I  <  t— — i —  ( I  v*  -  v  I  +  h  (  1 8  I  I  f  *  -  f  I  +  ,  . .  +  1 8   If*  -f  \))i 

ul       yll   -   1  -  |B   |hn  uy0  y0l        avlP0l  lxo   10l         |Pk'  '  -k   -kUh 

from  now  on.  we  suppose  that  h  <  h  =  -rr-r — -r—  .  let  B  =  6_  +  \b.\    +    ...    +    \&,  \ ; 

—     2 1 6      1 0  *  '  ^0 '    "i. '  k '  ' 

since  y*  =  y  _,  we  have 

|y*  -  yj  <   2hB| |u*  -  uQ| |    and     | f *  -  f±\    <  2hftB| |u*  -  uQ | |;    (k) 

in  the  same  way,  we  get  for  the  second  step: 

|y*  -  v  I  <  2fl v*  -  v   +  h( Ib   I  f  *  -  f   +  I B  I  If  *  -  f  I  +     +  I B  I  If*   -f   1)1 
u2       •y2l  -   u,yl   yl'    "^IPo111!   Ii"  +  -lpl»  1*0   I0I         IPkM  1-k    l-k" 

by  (h): 

|y*   -   y2l    <  2{2hB-||u*   -   u0||    +  h(2hBa||uJ   -  uQ|  |    +   b||u*   -   ujl)}; 
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\y*   -  y2|  <  Bh(3  +  2hn)|  |u*  -  uQ\\; 


since  B  and  ft  are  constant,  we  can  write 


y*  -  y0|  :  0(h) ||u* 


u.. 


2   j^    _   v-/i  !-0    0 


t 


repeating  the  same  process  k  times,  we  finally  get 

|y*  -  y.|  =  0(h)|  |u*  -  uQ|  |,      i  =  0,1,,..,*,  (5) 

|f  J  -  f.|  =  0(h)|  |ug  -  uQ|  |,       i  =  0,1,.  ..,k.  (6) 

Let  R(x)  and  R*(x)  be  respectively  the  polynomials  of  degree  k  such  that 

R(x.)  =  f .,      i  =  0,1,  ...,k, 


R*(x.)  =  f *,      i  -  0,1, ...,k; 


Pl'  '"'    Pk+2'  pl'  '••'  Pk+2'  comP°nents  respectively  of  k(uq)  and  k(u*)  are 
given  by 


px  =  pJ  =  y0, 


p2  =  R(xQ)  =  p*  =  R*(xQ)  =  fQ, 


P3  =  R(x_1);      p*  =  R*(x_3), 


Pk+2  ■  R(x   )j      PJ+2  =  R*(x_k); 


then,  there  exist  the  constant  coefficients  (independent  of  h)  such  that 
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p.  =  Z  a.  .f ., 

1   j=0  1J  J 


k 
p*  =     Z  a.  .f*: 

1   j=o  1J  J 


by  (6) 


p.  -  p*|  =  0(h)||u*  -  u0|.|,      i  =  0,1,. ..,k 


|R(u*)  -  k(uq)| I  =  0(h)| |ug  -  uQ| 


this  is  the  desired  result  since  for  h  small  enough  0(h)  <  1. 

Remark  that  we  have  proved  so  far  only  the  convergence  of  the 
polynomials  Pn  .  (x);  "but  the  same  arguments  apply  also  for  P,  .  (x)„ 


Proof  of  Theorem  1 

We  have  already  proved  the  convergence  of  the  simplified  starting 
procedure  and  of  the  starting  procedure. 

Let  us  call  U(x)  the  limit  polynomial  for  the  simplified  starting 
procedure;  by  construction  U(x)  satisfies  (2)  and  (3)  for  i  =  0;  it  remains 
to  prove  (3)  for  i  =  1,2, ...,k.   Let  y  ,  y  ,  ...,  y  be  the  values  obtained  by 

,  ,  12  K. 

integrating  with  the  k   Nordsieck  formula  and  starting  polynomial  U(x). 
Because  U(x)  is  the  limit  polynomial,  it  follows  that 

U'(x.)  =  f(x.,y.),     i  =  0,l,...,k;  (7) 

since  the  k   Nordsieck  formula  is  equivalent  to  the  k   Adams -Moulton  formula, 
it  follows  that  y  =  R(x  ),  where  R(x)  is  the  polynomial  of  the  degree  k+2 
satisfying  the  relations 

R(*0)  =  y0  =  u(x0^ 
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R'(x.)  =  U'(x.),      0,-1,. ..,-k, 


R'(x1)  =  f(Xl,yi); 


^y  (T);  R'OO  =  U'(x  ),  and  this  implies  that  R(x)  =  U(x),  y  =  U(x  ),  which 
proves  (3)  for  1=1,   For  i  =  2,3* .  ..>k  the  proof  follows  the  same  pattern. 

We  now  discuss  briefly  the  part  of  theorem  1  which  concerns  the 
starting  procedure.   Let  Pn(x)  and  Qn(x)  be  the  limit  of  the  polynomial  P   .(x) 
and  Q   (x).   We  know  that  Pn(x)  and  Qn(x)  are  independent  of  P  n(x);  let  us 
choose  P  n(x)  =  U(x),  where  U(x)  is  the  polynomial  which  satisfies  the  relations 
(2)  and  (3)  (its  existence  has  been  established  above.)/  by  arguments  similar  to 
those  used  previously,  it  follows  easily  that  Q_  n(x)  =  U(x)  =  Pn  n(x)  and 
consequently  Pq(x)  =  ^n^x^  =  u(x)- 

In  order  to  avoid  fastidious  repetitions,  we  renounce  to  give  the 
details  of  the  proof  of  theorem  2. 


Starting  Errors  in  the  Nordsieck  Method 

From  theorem  1,  it  follows  that  the  limit  initial  polynomial  U(x)  is 
the  same  for  the  starting  and  the  modified  starting  procedure.   Let  P  (x), 
P  (x),  ...  be  the  polynomials  corresponding  to  the  steps  0,1,2,.../  we  know  that 
U(x)  =  P0(x)  =  Pk(x)  (and,  in  fact,  also  U(x)  =  P  (x ),  P  (x),  .  .  .  ,P  ,(x)).   Let 
z(x)  be  the  exact  solution  of  the  differential  equation  with  z(x  )  =  y  ;  we  want 
to  apply  the  asympototical  theorem  of  Section  2  with  q  =  k/  we  have  to  determine 

U(xk)  -  z(xk),  (8) 


U'(x.)  -  z'(x.),      i  =  0,'l,...,k.  (9) 

Let  Q(x)  be  the  polynomial  of  degree  k+1  (as  U(x))  satisfying  the  relations 

Q(x0)  =  y0,  (10) 


Q'(x.)  =  z'(x.),      i  =  0,1,.  ..,k;  (11) 


■kQ- 


let 


R.  =  z(x.)  -  Q(x.),      i  =  0,1,  ...,k;  (12) 


we  have  the  relations 


U(xQ)  =  y 


U'(x.)  =  f(x.,U(x.)),      i  =  0,1,  ...,k,  (13) 


Q(xQ)  =  y0, 


»'  (x.)  =  f(x.,Q(x.)  +  R.),      i  =  0,1,  ...,k;         (lV 


let 


H(x)  =  U(x)  =  Q(x),  (15) 


t.  -  H'(x.),      i  =  0,1,  ...,k;  (16) 


since  H(xn)  =  0,  there  exist  the  constants  a.,  independent  of  h  such  that 

k 
H(x)  =  h  Za  t,     1  =  0,1,...,^  (17) 

j=0  1J  J 


we  define  the  quantities  g.,  i  =  0,1, ...,k  by  the  relations 

f(x.,U(x.))  -  f(x.,Q(x.)  +  R.)  =  g.(u(x.)  -  Q(x.)  -  R.),      i  =  0,1,.  ..,k; 

(18) 

if  U(x. )  -  Q(x. )  -  R.  =  0,  we  set  g.  =  0, . 

i     ^  i7    i    *  &i 

We  know  that  |g.|  <  Q;  subtracting  (lh)   from  (13),  we  get  by  (15), 
(16),  (17),  (18), 
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k 

t.  =  g.h  2  a  tj  -  g.R.,      i  =  0,1,  ...,k;  (19) 

J  —  \J 


this  is  a  linear  system  for  the  t.'s  if  we  consider  R.  as  known;  we  can  regard 
t.  as  function  of  h;  suppose  that  there  exists  a  constant  m  such  that 


R.  =  0(h  ),  i  =  0,1, ...,k;  it  follows  the  relati 


on 


1 


t.  =  -g.R.  +  0(hm+1).  (20) 

111  v  ' 


We  only  sketch  the  proof  of  (20 ).   Since  g.  is  bounded,  a.  .  independent  of  h, 

(20)  follows  from  (19),  if  we  show  that  t.  =  0(hm),  i  =  0,1, ...,m.   Let  t  be  the 

-> 
vector  of  components  t.,  w  the  vector  of  components  -g.R.  and  B  the  matrix  of 


elements  g.a..;  (19)  becomes 
1  ij7 


(I  -  hB)t  =  w;  (21) 


we  introduce  the  following  notations:   for  any  matrix  C,  the  norm  | |c| |  of  C 
is  the  sum  of  the  absolute  value  of  its  elements;  norms  of  vectors  and  matrices 
have  the  properties: 

a)  I |Cx| I  <  I |c| I  I |x| |; 

b)  I  |c  +  D|  I  <  I  |c|  I  +  I  |d|  I; 

c)  ||cd||  <  ||c||  •  ||d||; 

by  the  Gerschgorin  theorem,  if  h  <   1 11 1  ,  all  the  eigenvalues  of  hB  lie  in  the 

—  2  I  I  B  \   I 

circle  of  radius  0.5  around  the  origin;  the  matrix  I  -  hB  is  regular  and  its 
inverse  possesses  the  convergent  development 


(I  -  hB)"1  =  EhV; 
n=0 


by  b )  and  c ) 


00 


n=0 
-50 


by  (21)  and  a) 


—  x  -  n  ±5         — 


Now  we  can  determine  the  errors  (8)  and  (9);  "by  (12)  and  (15),  we  have 


U(x.)  -  z(x.)  =  H(x.)  -  R.,  i  =  0,1,... ,k; 


by  (17)  and  (20 ) 


m+2 


U(x  )  -  z(x  )  =  -h  Ea  g  R  -  R  +  0(h    ),      i  =  0,1,  ...,k;   (22) 


by  (ll),  (16)  and  (20 ) 


U'(x.)  -  z'Oc.)  =  t  =  -g.R.  +  0(hm+1),      i  =  0,1,  ...,k.      (23) 


In  order  to  get  more  useful  expressions,  let  us  suppose  (as  in  the  asymptotical 

df 
theorem)  that  G(x,y)  =  —  exists  and  satisfies  the  following  requirements: 

l)   |G(x,y*)  -  G(x,y)|  <  ty|y*  -   y|,      \|/  constant,  b  <  x  <  c; 


2)   G(x,z(x))  possesses  a  continuous  derivative j 


it  follows  that  g.  =  G(x.,y*),  where  y*€[U(x.),  z(x.)];  by  (22),  U(x, 
0(h  ),•  because  of  the  above  two  requirements,  we  have 


z(x.) 


g.  -  G(x.,z(x.))  +  0(hm)  -  G(xQ,z(x0))  +  0(h)  +  0(hm) 


Replacing  in  (22)  and  (23),  we  get  the  final  formulas 


*\ 


U(x.)  -  z(x.)  = 


U'(x.)  -  z'(x.)  - 

1        1 


■hG(xn,z(xn))  Ea  R   -  R  +  0(hm+2), 
j=0  J  J 


0'  v  0' 


>    i.  -   0,1,  .  .  .,k 


■G(x0,z(x0))R.  +  0(hm+1); 


J 


(2k) 
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We  recall  that 


G(x,y)  -  21^1 


Q(x)  is  the  polynomial  of  degree  k+1  which  satisfies  the  relations 


Q(xQ)  =  zQ,      Q'(x.)  =  z'(x.),      i  =  0,1,..  .,k; 


R.  =  z(x.)  -  Q(x.); 


a. .  are  the  coefficients  of  the  relations 


K(x  )  -  K(x  )  +  h  Za  K'(x  ),   i  =  0,l,  ...,k,  (25) 

j=0  J 


where  K(x)  is  any  polynomial  of  degree  k+1, 


Remark 


We  have  considered  the  errors  of  the  starting  procedure  for  the  limit 
L  U(x); 
some  j;  we  write 


polynomial  U(x);  but  practically,  we  must  deal  with  the  polynomial  P   .  for 


P.   (x.  )  -  z(x  )  =  (P   (x  )  -  U(x  ))  +  (U(x  )  -  z(x  )),      i  =  0,1,.  ..,k; 

by  choosing  P  n(x)  =  yn  +  fn(x  -  x„),  one  gets  easily  from  the  proof  of  the 
convergence  that 


P   (x  )  -  U(x  )  =  0(h2+J),      i  =  0,1,  ...,k.  (26) 

This  result  is  valid  for  both  the  starting  procedure  and  the  modified  starting 
procedure;  for  the  simplified  starting  procedure,  there  is  no  reason  for  thinking 
that  the  order  of  the  error  in  (26)  may  be  greater  than  h   ,°  however,  it  is 
very  likely  that  for  the  starting  procedure  we  have 
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P   (x  )  -  U(x  )  -  0(h2+2j),   i  =  0,1,  ...,k; 


no  proof  of  this  statement  has  "been  established;  if  that  is  true,  it  means 
that  the  speed  of  convergence  in  comparison  with  the  number  of  steps  of 
integration  to  execute  is  the  same  for  the  starting  and  the  modified  starting 
procedure;  the  simplified  starting  procedure  has  the  advantage  of  using  only 
forward  integration  which  remains  closer  to  the  true  solution  than  the  backward 
integration. 

Example 

We  consider  the  starting  procedure  for  the  first  Nordsieck  method  and 


suppose 


1)  enough  iterations  have  been  executed  for  using  (2k) }    when  U(x) 

is  replaced  by  P   .(x), 
V>  J 

2)  all  the  assumptions  necessary  to  the  application  of  the 
asymptotic  theorem  are  filled.   A  little  arithmetic  leads  to 
the  relations 


A  -  1 

a00  "  a01  "  °;   ai0  "  ail  "  2  ; 


R0  =  0;  Rx  =  -  X2  z'"(*0)  +  0(h^.); 


h3         „,,h 


y 


P0>.(Xl)-Z(xi)=^Z-(K0)+0(h') 


pi,j(xo'  ■  z(xo'  ■  °- 


if  u(x)  denotes  the  approximative  solution  obtained  with  this 
starting  procedure  and  the  first  Nordsieck  formula  and 
d(x)  =  u(x)  -  z(x),  we  have  by  the  asymptotical  theorem 
(with  q  =  1 ) 

d(x)  =  h3(ei(x)  +  0(h))  +  h3(eI];(x)  +  0(h)), 
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where 


§-  e_(x)  -  G(x,z(x))eT(x)  +  5J  z^\x),      e_(b)  =  0; 


dx   I  '        I        24 


^  en(x)  =  0(x,«(x))e  (x),      a  (b)  -  j|  z(3)(b); 


the  term  eJT(x)  represent  the  initial  error  due  to  the  starting  procedure, 


Errors  for  the  Starting  Procedures  for  the  Modified  Nordsieck  Formula 

All  we  have  said  above  can  be  repeated  word  for  word  for  the  k   modified 
Nordsieck  formula  except  that  we  must  replace  everywhere  k  by  k+1. 

h 

Example 

For  the  first  modified  Nordsieck  formula,  we  have: 

a00  =  a01  =  a02  =  °; 


_5         _8         zl 
al0  "'  12'   ail  "  12'   ai3  "  12' 


14         1 

a20   3'   a21  ""  3'  a22  "  y 


4  5 

R0  -  0,   Rx  =  1^  z(^(b)  +  0(h5),   R2  =  -  ^  z(5)(b)  +  0(h6); 


P0^.(x2)  -  z(x2)  =  -  -^  h5G(b,z(b))z(1+)(b)  +  ^  z(5)(b)  +  0(h6); 

p6,j(xo)  -  z,(xo}  =  °'        po,j(xi}  -  z'(xi}  =  -G(b>z(b))  h  z(i+)(b)  +  0(h5)' 

if  u(x)  denotes  the  approximative  solution  and  d(x)  =  u(x)  -  z(x): 
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d(x)  =  h3(ei(x)  +  0(h))  +  h5(ei;[(x)  +  eTI1(x)   +  0(h)), 


where 


|j  ex(x)  -  G(x,z(x))e].(x)  +  g|  z(J+)(x),      ej(h)  -  0; 
^  eIX(x)  =  G(x,z(x))eiI(x),      en(b)  =  -  j|  G(b,z(b))z(i+)(b)  +  ^  z(5)(b); 
^  eIlx(x)   =   G(x,z(x))eni(x),    em(b)  =  2B8  G(b,  z(b)  )z(U)(b ); 


the  starting  error,  represented  by  the  term  e   +  e    is  negligible  in  comparison 
with  the  term  eT« 
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SECTION  8.   CONCLUSION 

The  Nordsieck  method,  the  modified  Nordsieck  method  and  all  Nordsieck- 
like  formulas  considered  in  Section  k   look  very  elegant. 

Let  us  consider  once  more  the  matrix  formulation  of  the  multistep 
method  (3),  (k),    (5)  and  its  Nordsieck-like  equivalent  (6),  (7),  (8).   The 
classical  multistep  method  is  such  that  A  and  a  in  (k)   are  the  simplest  possible, 
i.e.,  each  normal  step  involves  the  least  possible  arithmetics.   If  v  is  the 
vector  of  the  iterative  process  (7)  and  v*  the  vector  obtained  after  the  modifi- 
cation of  h,  the  Nordsieck-like  method  is  such  that  the  matrix  D  which  transforms 

— »■     — * 

v  in  v*  is  diagonal;  the  corresponding  matrix  for  the  classical  multistep 

method  has  not  this  property. 

In  view  of  the  fact  that  keeping  the  step  constant  is  the  rule  and 
modifying  it  is  the  exception,  one  may  prefer  the  old  and  classical  formulation 
of  the  multistep  method. 
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SECTION  9-   TABLE  OF  RESULTS 

For  the  Nordsieck  method  and  the  modified  Nordsieck  method  (k  =  1,2, 3**0 
we  give  successively: 

1)  The  Nordsieck  formulation  of  the  method. 

2)  The  corresponding  multistep  corrector  formula: 

y(x  +  h)  =  y(x)  +  h(p_1f(x  +  h)  +  ...  +  3kf(x  -  kh)); 
the  truncation  error  is  denoted  by 

T  =  ChP+12(p+l), 

meaning  that  if  z(x)  is  a  function  which  possesses  (p+l)  continuous 
derivatives,  then 

z(x  +  h)  =  z(x)  +  h(p  z'(x  +  h)  +  ...  +  az'(x  -  kh))  +  T, 

T  =  ChP+1z(p+l)(|  ),    x  -  kh  <  I  <  x  +  h. 

3)  The  corresponding  multistep  predictor  formula  with  the  truncation 
error. 

k)     The  error  involved  in  the  modification  of  the  length  of  the  step. 
Suppose  the  differential  equation  y'  =  f(x,y)  is  integrated  by  the 
method  we  are  considering  with  the  step  h*  for  x  <  b  and  the  step 
h  for  x  >  b.   Let  u(x)  be  the  approximative  solution,  z(x)  be  the 
exact  solution  such  that  u(b)  =  z(b)  and  d(x)  =  u(x)  -  z(x)  for 
x  >  b;  the  error  satisfies  the  asymptotic  relation 

d(x)  -  hP(s(x)  +  0(h))  +  hq(t(x)  +  0(h)); 

the  term  s(x)  represents  the  accumulation  of  the  truncation  errors; 
the  term  t(x)  is  due  to  the  modification  of  the  step;  in  the  defini- 
tions of  s(x)  and  t(x)  we  shall  use  the  notations 
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(x)  .  M^l/y  =  z(x).  (1) 


h* 

p  =T' 


5)  The  error  involved  in  the  starting  procedure.   The  differential 

equation  y1  =  f(x,y)  is  integrated  for  x  >  b  with  prescribed  initial 
value  at  x  =  b;  if  u(x)  is  the  approximative  solution,  z(x)  the 
exact  solution  such  that  u(b)  =  z(b),  d(x)  =f   u(x)  -  z(x);  then  for 
x  >  b 

d(x)  -  hP(s(x)  +  0(h))  +  hr(w(x)  +  0(h)); 

the  first  term  represents  the  accumulation  of  the  truncation  errors 
and  the  second  the  error  due  to  the  starting  procedure.   The 
notation  (l)  will  be  used. 


Nordsieck  Formula  k  =  1 


y(x  +  h)  =  y(x)  +  h(f (x)  +  a(x)  +  -|  (f(x  +  h)  -  fP)}, 


fP  =  f(x)  +  2a(x), 


.(x  +  h)  =  a(x)  +  I  (f(x  +  h)  -  fP. 


Corrector  Formula 


y(x  +  h)  =  3I  (5f(x  +  h)  +  8f(x)  -  f(x  -  h)); 


1  ~2%   z 
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Predictor  Formula 


y(x  +■*)•  y(x)   +  -   (3f(x)   -   f(x   -   h)); 


1  i2   z         • 


Modification  of  h 


d(x)  =  h3(s(x)   +  0(h))  +  h3(t(x)  +  0(h)); 
s'(x)   -  g(x)s(x)   +  -^  Jk)(x),  s(b)   =   0; 

f(x)   =»  g(x)t(x),  t(b)   -i^£  z(3)(b). 

Starting  Procedure 

d(x)  -  h3(s(x)  +  0(h))  +  h3(w(x)  +  0(h)); 

s'(x)  =  g(x)s(x)  ±  ^  z{k)(x),  s(h)   =   0; 


w'(x)  =  g(x)w(x),      v(b)  =~   z(3)(b). 


Modified  Nordsieck  Formula  k  =  1 


y(x  +  h)  =  y(x)  +  h{f(x)  +  a(x)  +  b(x)  +  -|  (f(x  +  h)  -  fP)}, 


fP  =  f(x)  +■  2a(x)  +  3b(x), 

a(x  +  h)  =  a(x)  +  3b (x)  +  ^  (f(x  +  h)  -  fP), 

b(x  +  h)  =  b(x)  +  |  (f(x  +  h)  -  fP). 
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Corrector  Formula 


y(x  +  h)   =  j|  (5f(x  +  h)  +  8f(x)   -   f(x   -  h)); 


—  £-w- 


Predictor  Formula 


y(x  +  h)  =  -|  (23f(x)   -  l6f(x  -  h)  +  5f(x  -  2h)); 


,.#.»>. 


Modification  of  h 


d(x)  =  h3(s(x)  +  0(h))  +  h  (t(x)  +  0(h)); 


(x)  =  g(x)s(x)  +^z(^(x),  s(b)  =  0; 


t'd)   -  g(x)t(x),  t(b)   «   -<2p2   -Tf  +  H  .<«(»), 


Starting  Procedure 

d(x)  =  h3(s(x)  +  0(h))  +  h5(w(x)  +  0(h)); 


(x)   =  g(x)s(x)  +-i  z(U)(x),  s(b)  =  0; 


2K 


w'(x)   =  g(x)w(x),  w(b)  =  g  g(b)z(4)(b)  +  ^  z(5)(b), 
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Nordsieck  Formula  k  =  2 


y(x  +  h)  =  y(x)   +  h{f(x)   +  a(x)   +  b(x)   +  |   (f(x  +  h)   -   fP)}, 

fP  =  f(x)   +  2a(x)   +  3b(x), 
a(x  +  h)   =  a(x)  +  3b (x)  +  |   (f(x  +  h)   -   fp), 
b(x  +  h)   =  b(x)   +  i   (f(x  +  h)   -   fP)„ 


Corrector  Formula 


h  , 


y(x  +  h)  =  y(x)  +  ^  (9f(x  +  h)  +  19f(x)  -  5f(x  -  h)  +  f (x  -  2b)), 

720  z 


Predictor  Formula 


y(x  +  h)  =  y(x)  +  -|  (23f(x)  -  l6f(x  -  h)  +  5f(x  -  2h)); 


T  =  ^z(^) 


Modification  of  h 


d(x)  =  h\s(x)  +  0(h))  +  h\t(x)  +  0(h)),- 


'(x)  =  g(x)s(x)  +  -|2z(5)(x),      s(b)  =  0. 


t'(x)  =  g(x)t(x),      t(b)  =^^  z(4)(b). 
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Starting  Procedure 

d(x)  =  h*(s(x)  +  0(h))  +  h5(w(x)  +  0(h)); 

s'(x)  =  g(x)s(x)  +  ^|§  z(5)(x),      s(b)  =  0; 
v'(x)  -  g(x)w(x),      w(b)  =  -  ^  g(h)z(4)(b)  +  ^  z(5)(b)„ 

Modified  Nordsieck  Formula  k  =  2 

y(x  +  h)  =  y(x)  +  h{f(x)  +  a(x)  +  b(x)  +  c(x)  +  |  (f(x  +  h)  -  fp)}, 

fP  =  f(x)  +  2a(x)  +  3b  (x)  +  kc(x), 


a(x  +  h)  =  a(x)  +  3b (x)  +  6c(x)  +  j|  (f(x  +  h)  -  fp), 


b(x  +  h)  =  b(x)  +  i+c(x)  +  |  (f(x  +  h)  -  fp), 


c(x  +  h)  =  c(x)  +  ^  (f(x  +  h)  -  fP) 


Corrector  Formula 


y(x  +  h)  =  y(x)  +  ^  (9f(x  +  h)  +  19f(x)  -  5f(x  -  h)  +  f(x  -  2h)). 


I9hf  (5) 
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Predictor  Formula 


y(x  +  h)  =  y(x)  +  ^  (55f(x)  -  59?  (x  -  h)  +  37f(x  -  2h)  -  9f(x  -  3h)); 

_  2$lh^  (5) 
1  '"  720  Z    ' 
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Modification  of  h 

d(x)   =  h   (s(x)   +  0(h))   +  h5(t(x)   +  0(h)); 


(x)    =   g(x)s(x)   +  -&  z(5)(x),  s(b)   =   0; 


f(x)   =  g(x)t(x),  t(b)   =   -  £-    -   f2   +  1   z(5)(h), 


Starting  Procedure 

d(x)  =  h  (s(x)  +  0(h))  +  h5(w(x)  +  0(h)); 


(x)  =  g(x)s(x)  +  -g  z(5)(x),      s(b)  =  0; 


w'(x)  =  g(x)w(x),      w(b)  =  ^  z(5)(b). 


Nordsieck  Formula  k  =  3 

y(x  +  h)  =  y(x)  +  h(f(x)  +  a(x)  +  b(x)  +  c(x)  +  ||i  (f(x  +  h)  -  fP)} 

fP  =  f(x)  +  2a(x)  +  3b(x)  +  kc(x), 


a(x  +  h)  =  a(x)  +  3b(x)  +  6c(x)  +  i|  (f(x  +  h)  -  fP), 


b(x  +  h)  =  b(x)  +  kc(x)   +  i  (f(x  +  h)  -  fP), 


(x  +  h)  =  c(x)  +  -i  (f(x  +  h)  -  fP) 


2¥ 
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Corrector  Formula 

y(x  +  h)  =  y(x)  +  =|^(251f(x  +  h)  +  6k6f(x)   -  26kf(x  -   h)  +  106f  (x  -  2h)  -  19f  (x  -  3h)); 

1    ISo    ' 


Predictor  Formula 

y(x  +  h)  =  y(x)  +  g|  (55f(x)  -  59f(x  -  h)  +  37f(x  -  2h)  -  9f(x  -  3h)); 

_  2$lh^  (5) 
1    720  z    * 


Modification  of  h 

d(x)  =  h5(s(x)  +  0(h))  +  h5(t(x)  +  0(h)); 

s'(x)  =  g(x)s(x)  +  jgg  z(6)(x),      s(b)  =  0; 


f(x)  =  g(x)t(x),      t(b)  =  -  10P  j^g  -   9  z(5)(b) 


Starting  Procedure 

d(x)  =  h5(s(x)  +  0(h))  +  h5(w(x)  +  0(h)); 

s'(x)  =  g(x)s(x)  4-  ^  z(6)(x),      s(b)  =  0; 
V(x)  =  g(x)w(x),      w(b)  =  g§  z(5)(b). 
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Modified  Nordsieck  Formula  k  =  3 


251 


y(x  +  h)  =  y(x)  +  h{f(x)  +  a(x)  +  b(x)  +  c(x)  +  d(x)  +  |g  (f(x  +  h)  -  fP)}, 
fP  =  f(x)  +  2a(x)  +  3b  (x)  +  kc(x)   +  5d(x) 


,(x  +  h)  =  a(x)  +  3b (x)  +  6c(x)  +  10d(x)  +  ||  (f(x  +  h)  -  fP), 


b(x  +  h)  =  b(x)  +  kc(x)  +   10d(x)  +  ||  (f(x  +  h)  -  fp), 


c(x  +  h)  =  c(x)  +  5d(x)  +  5!  (f (x  +  h)  -  fP), 


d(x  +  h)  =  d(x)  +  ^  (f(x  +  h)  -  fp) 


Corrector  Formula 

y(x  +  h)  =  y(x)  +  =|rr  (251f(x+h)  +  6k6f(x)   -  26kf(x-   h)  +  106f  (x  -  2h)  -  19f(x-3h)). 


720 


T  -  -  4-   z(6) 
1     I£o  z 


Predictor  Formula 

y(x+h)  =  y(x)  +  T=!pr(l901f(x)  -  277^f(x-h)  +  26l6f  (x  -  2h)  -  127^f (x  -  3h)  ) 


720 

+  251f (x  -  kh)} 


T  -22L  (6) 
1  ""  288  z 


Modification  of  h 

d(x)  =  h5(s(x)  +  0(h))+  h6(t(x)  +  0(h)); 

-65- 


s'(x)    =   g(x)s(x)    +  rjk   z(6)(x),  s(b)    =    0; 


I5o 


4.1/    >         f    ur    \              4-^^          JJ8P11  -  7P2  -  180P  +  67     (6),,x 
t'(x)  =  g(x)t(x),  t(b)  = ^ zv    '(b), 


Starting  Procedure 

d(x)   =  h5(s(x)  +  0(h))   +  h7(v(x)   +  0(h)); 


(x)   =  g(x)s(x)  +  J-  z(6)(x),  s(b)   =   0; 


160 


w'(x)   =  g(xMx),  w(b)   =   -  j^_  g(b)z(6)(b)   +  J^  z(7)(b) 


Nordsieck  Formula  k  =  h 

y(x  +  h)   =   y(x)   +  h{f(x)   +  a(x)   +  b(x)   +  c(x)   +  d(x)   +  ^||   (f(x  +  h)   -   fP)}, 

fP  =  f(x)  +  2a(x)   +  3b(x)  +  kc(x)  +  5d(x), 


a(x  +  h)  =  a(x)  +   3b (x)   +  6c (x)   +  10d(x)   +  |?   (f(x  +  h)   -   fP), 


b(x  +  h)  =  b(x)  +  4c(x)   +  10d(x)   +  ||   (f(x  +  h)   -   fP), 


c(x  +  h)   =   c(x)   +  5d(x)  +  jj|   (f(x  +  h)   -   fP), 


d(x  +  h)  =  d(x)  +  ^   (f(x  +  h)   -   fP) 
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Corrector  Formula 

y(x  +  h)  =  y(x)  +  j^   {^75f(x  +  h)  +  l427f(x)  -  798f(x  -  h)  +  482f(x  -  2h ) 

-  173f(x  -  3h)  +  27f(x  -  hh)}, 

_      863hl  (7) 
1   "  ~  60480  z   • 

Predictor  Formula 

y(x  +  h)  =  y(x)  +  ^|q  (I901f(x)  -  277^f(x  -  h)  +  26l6f(x  -  2h) 

-  127^f(x  -  3h)  +  251f(x  -  kh)) 


T_  95h^z(6) 
1    288  z 


Modification  of  h 

d(x)  =  h6(s(x)  +  0(h))  +  h6(t(x)  +  0(h)); 

S'(x)^(x)s(x)+»z(T)^      s(b)  =  0; 
f(x)  =  g(x)t(x),     t(b)  =  -  j^   (2i+0p^  -  35P2  -  1792P  +  I587)z(6)(h), 

Starting  Procedure 

d(x)  =  h6(s(x)  +  0(h))  +  h7(w(x)  +  0(h)); 

s'(x)  =  g(x)s(x)  +^g_z^)(x),      s(b)  =  Oo 
V(x)  =  g(x)w(x),      v(b)  =  -  J^  g(b)z(6)(b)  +  ^  z(7)(b)„ 
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Modified  Nordsieck  Formula  k  -  h 
y(x  +  h)  -  y(x)  +  h{f (x)  +  a(x)  +  b(x)  +  c(x)  +  d(x)  +  e(x)  +  ^||  (f(x  +  h)  -  fp)}, 
fP  =  f(x)  +  2a(x)  +  3b(x)  +  4c(x)  +  5d(x)  +  6e(x), 
a(x  +  h)  =  a(x)  +  3b (x)  +  6c(x)  +  10d(x)  +  15e(x)  +  ^  (f(x  +  h)  -  fp), 
b(x  +  h)  =  b(x)  +  4c(x)  +  10d(x)  +  20e(x)  +  |  (f(x  +  h)  -  fP), 
c(x  +  h)  =  c(x)  +  5d(x)  +  15e(x)  +  U   (f(x  +  h)  -  fP), 
d(x  +  h)  =  d(x)  +  6e(x)  +  ^  (f (x  +  h)  -  fp), 
e(x  +  h)  =  e(x)  +  -|-  (f(x  +  h)  -  fP)„ 


Corrector  Formula 

y(x  +  h)  =  y(x)  +  ^^  (U75f(x  +  h)  +  lk-27t(x)   -  798f(x  -  h)  +  482f(x  -  2h) 


173f(x  -  3h)  +  27f(x  -  kh)) 


_   863^1  Z(T) 
1  "  "  60480  z   ' 


Predictor  Formula 


y(x  +  h)  =  y(x)  +  ^q  (U227f(x)  -  7673f(x  -  h)  +  9^82f(x  -  2h) 
-  6798f(x  -  3h)  +  2627f(x  -  kh)   -  425f(x  -  5h)). 


_  19087h7  (7) 
1  "  60480  z 
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Modification  of  h 

d(x)   =   h6(s(x)   +  0(h))   +  h7(t(x)   +  0(h)); 
s'(x)   .  g(x)s(x)   t  ^g  z(7)(x),        s(b)   =   0; 

t'(x)  =  g(x)t(x),  t(b)   =   -  —jL  (i6p5   -   3P3   -   Mp  +  28)z(7)(b) 

Starting  Procedure 

d(x)   -  h   (s(x)   +  0(h))   +  h7(w(x)   +  0(h)); 

•7 

s'(x)   =  g(x)s(x)   +  fg^  z(7)(x),        s(h)   =   0; 


w'(x)   =  g(x)w(x),  w(b)   =I|g_  z^)(b)< 
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APPENDIX  1 


Let  x  be  a  vector  of  m  components  x, ,  x^,  ....  x  .  f(x)  be  a  vector 

e  1'   2'    '  m 

function  of  x  with  components  f  (x),  ...,    f  (x).   For  any  vector  zf 


m 

Z  I z . I  is  the  norm  of  z . 
1-1  X 


Theorem 

If  for  any  two  x  and  x*,  there  exists  a  constant  c  <  1  such  that 


||f(x*)  -  f(x)|  |  <  c||x*  -  x||,  (1) 

then 

a)  the  equation 

x  =  f(x)  (2) 

— * 
possess  one  and  only  one  solution  y; 

b)  for  any  xn,  the  sequence  x_,  x  ,    ...,    x  ,    ...  defined  by 

x  ,  =  f(x  )  (3) 

n+1      n'  v  ' 

converges  towards  y. 

Proof 

This  theorem  is  only  a  particular  formulation  of  classical  results 
in  Banach  spaces  theory.   The  proof  results  from  the  two  statements: 

a)  the  equation  (2)  has  at  most  one  solution; 

b)  the  sequence  defined  by  (3)  converges  towards  a  solution  of  (2). 

a)  Let  y  and  y*  be  two  solutions  of  (2).   By  (l): 

My*  -  y||  ■  l|f(y*)  -  ?(y)|  |  <  c||y*  -  y|  I, 
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and  consequently  | | y*  -  y| |  -  0. 

.  -» 

b)   Let  x   -, ,  x  ......  x    be  the  components  of  x  .   In  order  to 

'  n,  1'      n,2'  n,m  n 

prove  the   convergence,    it   suffices   to  show  that  for  every  i, 
1  s  1.    . ..,   m,    the   sequence  x      .    is   a  Cauchy  sequence.      Let 
p  >  q  >  N: 

i  i  I  i -*  ~ *    i  i  l  I ~*  ~~ *        II  I  I-*  ~ *        I  I  I  I ~*  — *    l  l 

|Xp,i    -Xq,i'    "    l|Xp-XqM    -    l|Xq-Xp-lM    +    '  'Vl   "   Xp-2 '  '    +°°°+    "Vl"^1' 

(k) 

by  (1): 

||Xg   -  x1||    =    ||f(x1)   -   f(x0)||   <  cl|x1  -  x0||, 

||x3   -  x2||    =    ||f(x2)   -  Z&JW    <  c||x2   -  xjl    <  c2!!^   -  x0||, 


|x      .  -  x    I I    <  c    j  Jxn    -   x„J I; 

1    n+1       nn   -        '  '    1  0 '  " 


replacing  in  (h).    we  get 


p-1  oo  n 

n    c     M-*    -> 


x   .  -  x   .   <   x   -  x     Z  c  <   x_  -  xJ    Z  c   =  - — 
p, i    q,  l  —    1    O11        —  '  '  1    O'1    „      1-c'l    0 
*>  *'  n=q  n=N 


since  c  <  1  and  | |x  -  x  | |  is  constant,  for  given  e,  it  is  possible  to 

N  such  that   x   .  -  x 
p,i    c 

Let  y  =  lira  x  :  by  (3) 


find  N  such  that  llx   .  -x   .11  <  €,  p,q  >  N. 

p,i    q,i'   - 


n 


y  =  lim  x  .  =  lim  f  (x  ), 
J  n+1        v  n'- 

n-x»        n-*» 


by  (l),  f(x)  is  continuous  and  consequently  lim  f(x  )  =  y,  q„eod, 

n-*» 
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APPENDIX  2.   SOME  EQUIVALENT  "METHODS' 


Let  y'  =  f(x,y)  be  a  differential  equation  to  be  integrated  at 
x„.  x,  ...,x,  ...,x  =  x_  +  nh.   We  consider  h  as  a  fixed  number  and  denote 


'0>    '  V 


n    0 


by  y  the  approximative  solution  at  x  .   We  define  two  processes  for  computing 


Method  1 


For  n  >  0: 


-vp-> 

y  ,  =  p  v  +  yf     . , 
Jn+1   *     n      n+1' 


v  .,    =  Av     +f  _,  a , 
n+1     n    n+1  ' 


f  .  =  f  (x   _,y  ,  ); 

n+1      n+1'  n+1  ' 


(1) 
(2) 
(3) 


A  is  a  matrix  of  order  m  with  constant  elements;  p  and  a  are  constant  vectors; 

— >  — * 

y   is  a  scalar  constant;  v„  is  the  initial  vector,,   For  given  v_,  this  method 


allows  to  compute  successively  y  ,  y  , 


.,    provided  that  the  non-linear  system 
of  equations  (l),  (3)>  has  always  a  unique  solution  (we  shall  always  suppose 
that  it  is  the  case). 


Method  2 


For  n  >  0: 


y  .  =  s  ,  +  t  , . 

Jn+1    n+1    n+1' 


GO 


Sn+1  =  Vn  +  Vn-1  +  •"  +  ajVj  +  Vn  +  Vn-1  +  ""  +  ¥n-k  +  6fn+l 

(5) 
t    =  0  for  all  n,  except  for  n  =  w,  w+1,  ....,  w+p,  where  w  =  max(0,k~j); 
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j,  k,  p,  Q!_,  ...,0L.,    3n,  .  ..,3  ,  5  are  the  constants  which  characterize  the 

method  2.   Suppose  that  j  >  -1,  k  >  -1,  p  >  -1;  if  j  =  -1,  there  are  no  a-terms 

in  (5);  if  k  =  -1,  there  are  no  (3- terms  in  (5);  if  p  -  -1,  all  t-terms  are  zero 

and  the  method  (2)  becomes  the  classical  multistep  method;  furthermore  a./  0 

if  j  >  0,  Bk  /  0  if  k  >  0;  tv+v  tw+2,  ...,,  tw+p+]_,  flQ|  s_i;  ...,  s_.,  fQ,    f_±, 

f   ,  ....  f  ,  are  the  initial  values:  let  u_  denote  the  vector  of  dimension 
-2}  -k  '0 

p+j+k+3,  the  components  of  which  are  the  initial  values. 


Equivalence 

A  particular  method  1  and  a  particular  method  2  are  said  equivalent  if 
there  exist  two  rectangular  matrices  P  and  Q  with  constant  coefficients  (i„e.; 
independent  of  f(x, yj,  u_,  v  )   such  that  if  v  =  Pu  for  any  given  u  or 
u  =  Qvn  for  any  given  v.,  then  methods  1  and  2  give  the  same  values  for  y  } 
n  =  1,2,.-.. 


Theorem  1 

a)  For  any  method  1,  there  exists  an  equivalent  method  2. 

To )  For  any  method  2,  there  exists  an  equivalent  method  1. 

c)  Two  methods  of  type  2  are  equivalent  if  and  only  if  all  their 
characteristic  parameters  are  identical. 

Theorem  1  has  the  following  practical  meaning:   all  methods  1  are 
"essentially  equivalent"  to  a  classical  multistep  method.   These  methods  1, 
which  are  not  rigorously  equivalent  to  a  classical  multistep  method  (i.e.,  which 
are  equivalent  to  a  method  2  with  p  >  0),  do  not  seem  to  have  much  interest. 
Therefore  all  "interesting  methods  l"  can  be  discussed  in  terms  of  multistep 
formulas . 

The  proof  of  theorem  1  which  is  rather  painful  and  long,  will  not 
be  given  here. 

We  now  consider  two  other  methods,  method  III  (multistep  method)  and 
method  IV. 
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Method   III 

For  n  >  0: 


y     ,    ■  a.y     +  a.  y     ,+...+  a .  y     .+  Bnf     +  0,  f     .   +  ..'.+  0 .    'f      .       +  7  f 
Jn+1  CTn  lJn-l  J   n-j        K0  n       Kl  n-1  j+p  n-j-p 


n+lJ 


fn+1  "  f(xn+l^n+l)' 


i.  D,  a  ,  ....  a  ,  B~»  ....  B.   3  7  are  the  constants  which  characterize  the 

o}     r>       q>  >        j^  t-Q,  ,       J+P 


method;  j  >  0,  a  /  0;  p  >  0;  if  p  >  0,  then  B.+p  /  0;  yQ,  y_±,    .,.,  yn_. 
•J-P 


f^,  f  ,,  ....  f  .'   are  the  initial  values 
0'   -1'     '   -J-T 


Method  IV 

For  n  >  0: 

Zn+1  =  Vo,n  +  Vl,n  +   "'   +  Vj,n  "   6j+lgn  '    &j+2Sn-l   ""  "'    "   6j+Pgn-p+l  +  7W 

gn+1  =   f^n+l^n+l) 

z^       .    =  <x.z_       +  anzn        +   ...   +  a.z.       +   (8.   +  y)g       , 
0,n+l  0  0,n  1  l,n  J   J,n        v   j  /ton+l^ 

z.         .    =   z.    .        +  5.    .g     .,  i  =  1,2, . „ „,j; 

i,n+l  i-I, n  j-i°n+l'  '    '        *°' 

j>  P»  Qvn,  •  ••>  Ct.,  5  ,  ...,  6.,  5.  , .  6.   ,  7  are  the  constraints  which 
o>    *>>      q>  >      j>   o'       J    J+l   J+P 

characterize  the  method;   j  >  0,  a.   £   0,  p  >  0;  if  p  =  0,  there  are  no  6- terms 

J 


in  the  expression  of  z^±;    if  p  >  0,  5J+p  /  0;  zQ  Q,  . ...,  z^.  Q,  gQ,  g_^  . 


1  ^ 


g    are  the  initial  values. 

"*-    IT 


Theorem  2 


If  the  coefficients  B_,  ...,  B.    of  method  III  and  the  coefficients 

0'       J+P 

S~f  &.,...,  5.   of  method  IV  satisfy  the  relations 
0*  1'    '  J+p                 J 
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J 

P.  =  Zq  6.  .    -  8.  ...      i  =  0,l,...,p+j;  (6) 

Ki    _   r  j+i-r    j+l+i 


where  one  defines  6.    .  =  0  for  i  >  0,  and  if  the  relations 

j+p+i 

Sn-i  =  fn-i'      1  =  OAj.-.iP-I^  (7). 


p+i 

z.    =  y   .  +  E  6.  .   f   ,      i=  0,1, ...,j,  (8) 

i,n   Jn-i     _  J-i+r  n-r^  '    '        '°3  v 

3  r=0 


are  satisfied  for  n  =  0,  then  they  also  hold  for  n  =  1.  2.  .  .  .  and  z  =  y  „ 

3  J      3  n    n 

The  proof  of  theorem  2  by  induction  with  respect  to  n  is  straightfor- 
ward and  will  not  be  given  here. 


Remarks 

1)  For  given  p.,  (6)  is  a  linear  system  of  equations  for  8„,  .  .  . ,  8 .   ;  as  it  is 

easy  to  show  the  matrix  of  the  system  is  triangular;  all  diagonal  elements 

are  equal  to  a.;  by  hypothesis  a.  f  0  and  the  matrix  is  regular.   This 

J  J 

establishes  the  equivalence  of  methods  III  and  IV. 

2)  One  iteration  in  methods  III  or  IV  involves  practically  the  same  number  of 

operations.   However  method  III  needs  2j  +  p  +  2  "past  values"  y  3    y   ,..., 

y   .,  f  ,...,f   .    (supposing  p  >  0),  but  method  IV  needs  only  j  +  p  +  1 

"past  values,"  z^   ,  z.,   .....  g  ,...,g  ,   .   That  curious  property  would 
'  ,  r0,  n'  '  IfXx'  n'    'ton+l-p  *     *       J 

have  perhaps  been  useful  in  a  time  where  the  computers  had  only  small 
memories.   One  can  prove  that  for  two  different  sets  of  initial  values  in 
method  IV,  the  corresponding  sequences  of  values  z  are  also  different; 


this  is  not  the  case  for  method  III,  since  the  equations  (7)>  (8)  considered 

as  a  linear  system  for  y,...,y   .,  f,...,f   .    have  more  unknowns  than 

n'    3dn.-y      n3        '  n-j-p 

equations;  however,  the  corresponding  matrix  has  the  maximal  rank, 
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